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Abstract 


Many  low  level  visual  computation  problems  such  as  focus,  stereo,  optical  flow, 
etc.  can  be  formulated  as  problems  of  extracting  one  or  more  parameters  of  a  non¬ 
stationary  transformation  between  two  images.  Because  of  the  non-stationary  nature, 
finite-width  windows  are  widely  used  in  various  algorithms  to  extract  spatially  local 
information  from  images.  While  the  choice  of  window  width  has  a  very  profound 
impact  on  the  quality  of  results  of  algorithms,  there  has  been  no  quantitative  way 
to  measure  or  eliminate  the  negative  effects  of  finite- width  windows.  To  address  this 
problem.  We  introduce  two  sets  of  filters,  “moment”  filters  and  “hypergeometric” 
filters.  Due  to  their  recursive  properties,  these  filters  allow  the  effects  of  finite-width 
windows  and  foreshortening  to  be  explicitly  analyzed  and  eliminated. 

We  develop  one  paradigm  to  solve  general  one-parameter  extraction  problems  us¬ 
ing  moment  filters,  and  another  one  to  solve  general  multiple  parameter  extraction 
problem  using  hypergeometric  filters.  We  apply  these  paradigms  to  problems  of  focus 
and  stereo,  in  which  one  parameter  is  extracted  at  every  pixel  location,  and  optical 
flow,  in  which  two  parameters  are  extracted.  We  demonstrate  that  our  algorithms 
based  on  moment  filters  and  hypergoemetric  filters  achieve  much  higher  precision 
than  other  techniques. 

Keywords;  Focus,  Stereo,  Optical  Flow,  Gabor  Filter,  Moment  Filter,  Hypergeo¬ 
metric  Filter,  Low-level  Processing,  Computer  Vision,  Image  Processing. 


1  Introduction 


Finite- width  windows  are  widely  used  in  various  vision  algorithms  to  extract  spatially 
local  information  from  images.  While  the  choice  of  window  width  has  a  very  profound 
impact  on  the  quality  of  results  of  algorithms,  there  has  been  no  quantitative  way 
to  measure  or  eliminate  the  negative  effects  of  finite-width  windows.  We  propose 
two  sets  of  filters,  “moment”  filters  and  “hypergeometric”  filters.  Due  to  their  re¬ 
cursive  properties,  the  effects  of  finite-width  windows  can  be  explicitly  analyzed  and 
eliminated. 

Many  low  level  visual  computing  problems  such  as  focus,  stereo,  optical  flow,  etc. 
can  be  formulated  as  problems  of  extracting  one  or  more  parameters  of  a  “nonsta¬ 
tionary”  transformation  between  two  images.  By  nonstationary  we  mean  that  the 
parameters  are  position  dependent  in  the  image.  For  example,  we  can  formulate  a 
general  correspondence  problem  as  one  in  which  we  extract  disparity  as  one  parame¬ 
ter  in  a  position-dependent  phase-shifting  transformation.  The  nonstationary  nature 
of  the  parameters  forces  us  to  use  finite- width  filters  for  extracting  information  from 
images.  Despite  the  fact  that  most  algorithms  compute  the  parameters  by  either 
explicitly  convolving  images  with  the  finite-width  filters  (e.g.  phase-based  stereo) 
or  implicitly  combining  the  outputs  of  finite-width  filters  in  the  spatial  domain  (e.g. 
SSD  in  stereo),  the  effects  of  finite  width  are  seldom  analyzed.  Most  of  times,  the  in¬ 
formation  extracted  from  finite-width  filters  are  regarded  as  approximations  of  those 
from  infinite-width  filters.  But  such  approximations  make  the  computations  based 
on  the  extracted  information  much  less  accurate  or  even  totally  impossible. 

The  most  fundamental  difference  between  finite-  and  infinite-width  filters  in  the 
Fourier  domain  is  that  the  bandwidth  of  finite-width  filters  must  be  non-zero  while 
that  of  infinite-width  filters  is  zero.  Since  a  nonstationary  transform  in  the  Fourier 
domain  can  be  an  arbitrary  curve  (or  a  surface)  within  the  passband,  the  approx¬ 
imation  in  ignoring  the  effects  of  finite  width  is  equivalent  to  using  a  constant  to 
approximate  the  arbitrary  curve  in  the  passband.  The  moment  filters  and  hypergeo¬ 
metric  filters  approximate  the  curve  by  a  polynomial,  whose  order  can  be  arbitrarily 
high  so  that  the  approximation  can  be  as  accurate  as  a  specific  task  requires.  There¬ 
fore,  the  extraction  of  parameters  based  upon  the  polynomial  approximation  is  much 
more  accurate. 

We  apply  the  technique  of  moment  filters  to  the  problems  of  focus  and  stereo. 
Both  the  focus  problem  and  stereo  problem  are  examples  of  the  general  one-parameter 
extraction  problem.  Traditionally,  people  estimate  the  parameter  by  comparing  the 
outputs  of  convolutions  of  the  two  images  with  the  same  filter,  usually  a  Gabor  filter. 
Our  moment  filter  approach  provides  a  new  insight  into  this  traditional  approach  by 
showing  that  the  output  of  the  convolution  of  one  image  with  the  Gabor  filter  can 
be  represented  as  an  infinite  sum  of  the  output  of  convolutions  of  the  other  image 
with  moment  filters.  From  this  point  of  view,  the  traditional  approach  is  simply 
the  zeroth  order  truncation  of  the  infinite  sum,  which,  in  many  cases,  is  the  major 
source  of  inaccuracy.  By  using  higher-order  moment  filters,  we  are  able  to  reduce  the 
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inaccuracy. 

When  extracting  parameters  by  using  finite- width  filters,  people  usually  assume 
that  the  parameters  are  constant  within  the  width  of  filters.  But  another  aspect  of 
the  nonstationarity  is  that  the  parameters  actually  change  within  the  effective  width 
of  the  filters,  which  we  refer  to  as  the  shift  variance  problem.  Such  problem  is  referred 
to  as  foreshortening  in  stereo,  and  affine  matching  problem  in  2D  image  registration. 
Traditional  infinite-width  filters  are  unable  to  deal  with  the  situation  because  the 
infinite-width  Fourier  transform  doesn’t  converge.  The  recursive  properties  of  the 
moment  filters  allow  modeling  this  effect  as  a  modification  of  the  Gabor  filter  con¬ 
volution  by  convolutions  of  the  first-order  moment  filters.  Applying  this  technique 
to  the  focus  and  stereo  problem,  we  demonstrate  that  we  can  not  only  eliminate  the 
negative  effect  on  the  estimation  of  the  parameter  caused  by  the  foreshortening,  but 
also  estimate  the  the  degree  of  foreshortening  quantitatively.  Such  a  quantitative 
measurement  provides  a  new  approach  to  the  shift  variance  problem. 

Both  the  traditional  approaches  and  the  moment  filter  approach  suffer  from  the 
frequency  sampling  problem,  i.e.  there  is  no  way  that  the  frequency  domain  can 
be  sampled  completely  and  nonredundantly.  If  the  frequency  domain  is  sampled 
sparsely,  much  information  is  abandoned,  while  if  it  is  sampled  densely,  the  results 
from  different  frequency  bands  are  highly  correlated,  which  makes  the  merge  of  the 
results  difficult.  Extending  the  idea  of  the  moment  filters,  we  developed  another  set 
of  filters,  “hypergeometric  filters”.  The  advantage  of  the  hypergeometric  filters  is 
that  they  provide  a  complete  and  non-redundant  decomposition  of  the  local  signal, 
even  though  they  sample  the  frequency  domain  in  a  nonuniform  fashion.  Based  on 
hypergeometric  filters,  the  general  problem  of  extracting  parameters  of  nonstation¬ 
ary  transformation  can  be  formulated  as  a  multidimensional  minimization  problem. 
Like  the  moment  filters,  the  hypergeometric  filters  also  have  recursive  properties  in 
both  the  spatial  and  frequency  domain,  and  therefore,  are  also  capable  of  extracting 
parameters  with  high  precision  and  dealing  with  the  foreshortening  problem. 

We  apply  the  hypergeometric  filter  technique  to  the  optical  flow  problem,  which 
can  be  formulated  as  extracting  two  parameters  from  a  nonstationary  transformation 
between  two  images.  A  2D  conjugate  gradient  minimization  algorithm  is  employed  to 
compute  the  optical  flow  between  two  adjacent  frames.  Not  only  can  this  algorithm 
produce  very  precise  estimation  of  the  optical  flow,  but  also  a  meaningful  error  esti¬ 
mation,  which  represents  the  aperture  problem  as  covariance  matrices,  is  computed 
based  on  the  minimization.  Finally,  we  quantitatively  compare  the  performance  of 
the  hypergeometric  filter  approach  with  other  popular  techniques. 

2  Related  Research 

The  paradigm  of  extracting  parameters  by  convolving  images  with  filters  which  are 
local  in  both  the  spatial  and  frequency  domain,  has  been  successfully  used  in  many 
low-level  vision  tasks,  such  as,  motion  analysis,  stereo  matching,  texture  analysis, 
and  focus  measure  in  literature.  Adelson  and  Bergen  [2]  and  Heeger  [10]  modeled 
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motions  in  2D  image  space  as  orientations  in  3D  spatiotemporal  space  by  introducing 
3D  oriented  filters  to  measure  image  velocities.  More  recently,  Fleet  and  Jepson  [7] 
modeled  the  normal  velocity  as  a  function  of  local  phase  changes,  and  then  used 
Gabor  filters  to  measure  changes  of  phase  at  every  pixel  location.  For  the  stereo 
matching  problem,  Weng  [27],  Sanger  [23],  Fleet  et  al.  [8],  and  Langley  et  al.  [16] 
proposed  using  filters  to  extract  phase  information,  then  computed  disparities.  Jones 
and  Malik  [12],  on  the  other  hand,  applied  a  set  of  linear  spatial  filters  to  images  and 
used  responses  from  those  filters  as  matching  features.  Because  the  focus  difference 
between  two  images  can  be  represented  as  a  blurring  convolution,  many  researchers 
have  also  proposed  using  spatial  frequency  analysis  to  compute  the  focus  difference 
[29,  4,  26,  21].  The  spatial  frequency  approach  also  achieved  great  success  in  texture 
segmentation  [5,  11,  15]  and  shape  recovery  from  texture  [14,  17]. 

The  problem  with  finite-width  windows  mentioned  in  the  previous  section  appears 
in  every  application.  Realizing  finite-width  filters  cause  errors  and  instability,  some 
researchers  found  ways  to  deal  with  this  problems  in  literature.  One  way  is  to  assume 
that  there  are  one  or  more  dominating  frequency  components  in  images.  Then  around 
the  dominating  frequencies,  the  effects  of  finite  width  are  negligible.  Though  this 
approach  is  effective  when  the  assumption  is  valid,  ordinary  images  usually  don’t 
contain  any  dominating  components.  Another  way  proposed  in  literature  is  using  two 
techniques,  a  stability  criterion  for  detecting  the  situations  where  the  filter  output 
is  dominated  by  window  contamination  [6,  30]  and  instantaneous  frequency  as  the 
spatial  derivative  of  phase.  In  relating  to  our  approaches  proposed  in  this  report, 
both  techniques  are  actually  special  cases  of  the  moment  filter  approach.  Both  the 
stability  criterion  and  instantaneous  frequency  can  be  directly  computed  from  the  first 
order  moment  filter.  From  this  point  of  view,  the  moment  filter  approach  generalizes 
these  two  techniques. 

Another  common  practice  in  literature  in  using  finite-width  filters  is  that  low 
frequency  components  are  usually  abandoned  due  to  two  reasons.  One  is  that  the 
low  frequency  components  are  usually  contaminated  by  the  DC  component  because 
Gabor  filters  have  non-zero  DC  bias,  and  the  DC  component  is  usually  very  strong. 
The  other  reason  is  that  the  negative  effects  of  finite  width  is  much  severer  when 
the  frequency  is  low  than  that  when  the  frequency  is  high.  Unlike  Gabor  filters,  the 
hypergeometric  filter  contains  no  DC  component  as  long  as  the  peak  frequency  is  not 
zero.  And  since  we  eliminate  the  effects  of  finite-width  windows,  the  low  frequency 
components  can  also  be  used  to  compute  high  precision  results.  We  will  show  that  the 
hypergeometric  filter  approach  provides  a  canonical  way  of  sampling  the  frequency 
domain  so  that  all  information  is  gathered,  including  low  frequency  components  which 
actually  contain  surprisingly  rich  information. 

3  Moment  Filter  Approach 

In  this  section,  we  will  propose  the  “moment  filter”  decomposition  of  the  signal  and 
demonstrate  how  this  decomposition  can  be  effectively  applied  to  problems  of  focus 
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and  stereo.  When  we  apply  such  a  decomposition  to  problems  of  focus  and  stereo,  we 
show  that  the  effects  of  finite  width  can  be  represented  by  high  order  moment  filters 
though  Taylor’s  expansion  of  the  transformation.  Furthermore,  we  are  able  to  solve 
the  foreshortening  problem  in  both  cases  by  using  the  spatially  recursive  property  of 
moment  filters,  while  the  traditional  Fourier  analysis  fails  because  the  Fourier  integral 
doesn’t  converge. 


3.1  Theory  of  Moment  Filters 

3.1.1  Definitions  of  Moment  Filters 


We  define  the  ith  order  moment  filter  in  the  spatial  domain  as 


in  which 


1  -f 

mfx)  =  pi{x)  _  e  2<t2  e 
y/zTra 

(1) 

'  \  i  =  2„  +  l, 

(2) 

where  n  is  an  integer,  j  is  -y/^,  and  L{n,  a,  x)  is  the  generalized  Laguerre  polynomial 
[20],  with  L(0,  c,x)  =  1.  Figure  1  illustrates  some  of  the  moment  filters  in  the  spatial 
domain. 

Apparently  the  zero-order  moment  filter  is  the  Gabor  filter  G{x-,  fo)  with  the  peak 
frequency  fo: 

.  (3) 

v27r  cr 

In  the  Fourier  domain,  the  moment  filter  can  be  defined  as 


M,(f)  =  nmi(x)\  =  (/  -  (4) 


Therefore,  the  peak  frequency  of  the  fth  moment  filter  is 

Peak[M,(/)]  =  /o±— , 
c 


(5) 


which  moves  away  from  the  original  peak  frequency  as  i  increases.  Figure  2  illustrates 
the  pass  bands  of  moment  filters  as  i  increases,  assuming  fo  =  0.  Note  that  only  zero- 
order  moment  filter  has  one  single  peak,  all  others  have  two  peaks. 

Another  way  to  understand  the  moment  filters  is  that  the  profiles  of  the  filters,  i.e. 
mfx)  when  fo  —  0,  can  be  generated  by  differentiating  the  Gaussian  function  g{x): 

^i£;s{x)]  =  ufrm=^)] = (6) 

From  Eq.  4  and  Eq.  1,  we  then  have 

mix)  =  (7) 
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Figure  2:  The  Pass  Band  Moves  as  i  Increases 


3.1.2  Decomposition  by  Moment  Filters 

In  this  section,  we  are  considering  moment  filters  based  upon  one  Gabor  filter,  which 
has  the  peak  frequency  /q.  We  also  limit  our  discussion  to  one  specific  spatial  location 
xo  =  0  unless  we  advise  otherwise. 

Let 

k(x)  =  p'(x)e"",  (8) 

and  the  definition  of  an  inner  production  over  the  function  space  as 


<  g{x),h{x)  >= 


y/^a 


f+CO  ^2 

/  g{x)h*{x)e  2a^dx, 
J  —  OO 


where  *  denotes  complex  conjugate,  we  then  have 


(9) 


1.  The  functions  ki  in  Eq.  8  are  orthogonal  with  each  other  under  the  definition  of 
the  inner  production,  i.e. 


<  ki{x),kj{x)  >—  0,  whem  ^  j. 


(10) 


2.  For  any  function  ({x)  which  has  finite  energy,  we  can  decompose  the  function 
into  a  weighted  sum  of  ki.  In  other  words,  the  set  of  function  ki  is  a  complete 
and  orthogonal  basis  set  as 


Wi  =  <  ^{x),ki{x)  >, 

(11) 

(12) 
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The  proof  can  be  found  in  Appendix  A. 

Let  us  look  into  the  inner  production  in  Eq.  11, 
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This  equation  can  also  be  represented  in  the  Fourier  domain  as, 

1 

Wi  =  J'[^(a;)]  (8)^[pi(a;)-=-e  2-^] 

V27rcr 


=  /'■”  /-[{(xjK/o  - 
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In  other  words,  the  coefficients  Wi  from  convolutions  of  the  signal  with  the  mo¬ 
ment  filters  can  be  used  to  losslessly  reconstruct  the  original  signal.  Therefore,  the 
decomposition  of  the  original  signal  by  moment  filters  is  complete  and  non-redundant. 
Eq.  13  and  Eq.  14  show  the  computation  of  w,  in  the  spatial  and  the  Fourier  domains. 
We  will  use  these  two  equivalent  representations  in  our  later  analysis. 

3.1.3  Important  Properties  of  the  Moment  Filters 

We  here  list  the  important  properties  of  the  moment  filters.  The  proof  of  these 
properties  can  be  found  in  Appendix  B. 

•  They  are  recursive  in  the  Fourier  domain: 

(/-/o)M,(/)  =  M,+i(/).  (15) 


•  They  are  recursive  in  the  spatial  domain: 


xmi{x)  =  —  a'^mi+-i{x)), 


where  j  =  \/— 1. 


•  They  are  recursive  with  respect  to  the  differential  operation: 


j{mi+i{x)  -  fomi{x)). 
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•  The  instantaneous  frequency  is  defined  as  the  spatial  derivative  of  the  phase  in 
[23,  6].  Usually  replacing  the  peak  frequency  of  the  Gabor  filter  by  instanta¬ 
neous  frequency  will  increase  the  accuracy  as  claimed  in  [6].  The  instantaneous 
frequency  can  also  be  represented  by  moment  filters  as 

^«x„)  =  /„  -  Re  (^)  .  (18) 

where  Re  means  the  real  part,  ^  is  the  phase,  and  Wo  and  Wi  are  coefficients  in 
Eq.  13. 

•  The  stability  criterion  is  proposed  by  Fleet  et  ah  in  [6]  to  locate  where  the  infor¬ 
mation  extracted  by  finite-width  filters  is  substantially  different  from  those  by 
infinite-width  filters.  In  other  words,  the  extracted  information  is  overwhelm¬ 
ingly  contaminated  by  window  effects  when  the  stability  criterion  is  not  satisfied. 
We  can  represent  the  stability  constraint  as 


T=| 


1  da  ,,, 

-T-  '  + 

a  dxo 


r  ||2^  II  ^^1  IP 
dxo  °  II  Wo  IP 


(19) 


The  three  recursive  properties  establish  the  interdependence  among  moment  filters. 
They  are  the  most  important  ones.  The  last  two  properties  demonstrate  that  both 
instantaneous  frequency  and  stability  criterion  techniques  are  ad  hoc  utilization  of 
information  from  the  first-order  moment  filter. 


3.2  Moment  Filters  For  Focus 

3.2.1  Problem  Definition 


The  method  of  depth  from  defocus  is  to  obtain  depth  information  using  a  quantitative 
measure  of  difference  of  focus  between  two  images.  For  simplicity,  we  assume  a  Gaus¬ 
sian  model  of  the  blurring  function,  even  though,  as  we  will  see  later,  the  approach  is 
applicable  to  any  model.  Let  us  denote  two  images  as  fo(a^)  and  in  the  spatial 

domain  and  /o(/)  and  /i(/)  in  the  Fourier  domain.  The  relations  between  these  two 
images  are: 


Hf) 


1  f+oo  1 

io{x)  ®  —j= — e  =  /  «o(f)-^= — e  di 

y^l'KO'Q  J—oo  Y27r(7o 


Hf) 


-/“-S/a 


(20) 

(21) 


The  difference  of  focus  between  two  images  is  characterized  by  ao.  If  wq  is  not  a 
constant  within  the  window  but  changes  linearly  in  neighborhood  of  window  center 
xq  (referred  as  e(xo)),  i-e 

(To  =  5(1  -I-  pix),  X  e  t{xo),  (22) 


8 


we  can  expand  the  second  integral  term  in  Eq.  20  as  follow: 


1  1  {x  -  ty  -  5^ 

-  (23) 

where  we  truncated  terms  of  higher  order  of  fix  based  upon  the  assumption  that  [i  is 
small. 

Now  we  need  to  solve  the  following  two  problems: 

•  Finite  Window  Problem:  Since  Eq.  21  is  based  upon  an  infinite  window,  how 
can  we  compute  (Tq  when  the  window  is  actually  finite? 

•  Shift  Variance  Problem:  Since  Eq.  21  is  based  on  the  assumption  that  (Jq  is  a 

constant  within  the  window,  how  can  we  compute  s  in  Eq.  22  when  (Jq  is  not  a 
constant,  namely  //  0? 

Another  way  to  understand  Eq.  22  and  23  is  that  s  represents  the  blurring  differ¬ 
ence  at  the  pixel  and  ft  the  gradient  of  the  blurring  difference  at  the  pixel.  When 
the  blurring  difference  s  is  converted  into  depth,  the  gradient  fi  is  converted  into 
the  surface  gradient,  or  surface  tilt.  Therefore,  solving  the  shift  variance  problem  is 
equivalent  to  simultaneously  computing  depth  and  shape.  In  the  following  discussions 
we  will  only  use  a  single  Gabor  filter  with  the  peak  frequency  fo  and  the  moment 
filters  based  on  the  Gabor  filter.  In  practice,  however,  we  can  optimally  merge  results 
from  multiple  Gabor  filters  by  Kalman  filtering. 


3.2.2  Finite  Window  Problem 


In  this  section,  we  will  our  discussion  to  the  situation  when  ^  =  0,  i.e.  CTo  is  a 
constant  within  the  window.  Let’s  define  Ui  and  Vi  as  the  convolution  of  the  first  and 
the  second  image  with  the  zth  order  moment  filter,  i.e. 


U.  = 

Vi  = 


/+00 

iQ{x)mi{x)dx 

'OO 

/+00 

Io{f)Mi{f)df, 

■OO 

/  +  00 

ii{x)mi{x)dx 

-OO 


/+00 

h{f)Mi{f)df 

-OO 

(-!)■'  I*"  Io(f)e-r-t/^M,{f)df. 

J—OO 


(24) 


(25) 


By  Taylor’s  expansion,  in  the  neighborhood  of  /o,  i.e.  e(/o),  we  have, 

e-r-ln  «  ,  /  g  (26) 
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where  higher  order  expansions  are  truncated. 

Replacing  Eq.  26  into  Eq.  25,  and  using  the  recursive  property  of  moment  filters 
in  the  Fourier  domain,  we  obtain. 


^0  ^  J_  Hf)  -  fo(^Uf  -  /o) 


^0^(1  -  fyo) 


if  -  foY 


-/o"^oV2 


(27) 


Therefore,  we  have  the  equation  to  compute  cTq, 


^^(1  -  fl^l) 


al  =  ^  (^ln(f/o  +  -  -  2 


U2)-\nVo 


Compared  with  the  equation  from  Eq.  21,  i.e. 


(28) 


<^o  =  72(lnM/o)-lnA(/o)), 

JO 


(29) 


we  can  see  that  if  we  make  the  assumption  that  the  effects  caused  by  window  are 
negligible,  i.e. 

,  Co  loifo) 


(30) 


this  may  result  into  substantial  errors  if  either  Ui  or  U2  happens  to  be  not  very  small 
comparing  to  Co-  ■ 

A  more  intuitive  way  of  looking  at  the  difference  between  Eq.  28  and  Eq.  29  is  illus¬ 
trated  in  Figure  3.  The  solid  line  is  the  transformation  curve  e~-^  with  ao  =  1.0. 
The  peak  frequency  /o  =  0.8,  and  the  passband  of  the  Gabor  filter  is  (0.4, 0.8).  The 
dotted  curve  is  the  zeroth  order  approximation,  i.e.  Eq.  29  is  equivalent  to  using 
a  constant  to  approximate  the  transformation  curve  within  the  passband.  The  first 
order  approximation  (dash-dot  curve),  i.e.  Eq.  28  with  only  the  linear  term,  is  equiv¬ 
alent  to  using  a  line  to  approximate  the  transformation  curve  within  the  passband. 
The  second  order  approximation  (dashed  curve),  i.e.  Eq.  28,  is  equivalent  to  using  a 
quadratic  curve  to  approximate  the  transformation  curve  within  the  passband.  Gen¬ 
erally,  if  we  use  up  to  the  nth  order  moment  filter,  in  the  Fourier  domain,  we  are 
actually  fitting  the  transformation  curve  by  an  nth  order  polynomial.  Certainly,  it 
will  drastically  decrease  errors  comparing  to  the  constant  approximation. 

Now,  let  us  look  back  the  truncation  in  Eq.  26.  We  truncated  terms  of 
for  i  >  2.  Correspondingly,  we  truncated  terms  of  ^  for  f  >  2  in  Eq.  27.  These 
truncations  are  valid  only  if  ^  generally  decreases  in  magnitude  as  i  increases.  Since 
Ui  is  generated  by  convolving  images  with  mfx)  as  in  Eq.  13,  the  expected  magnitude 
of  Ui  will  be  proportional  to  the  square  root  of  the  energy  of  the  filter: 


Ei  =  II  mfx)  \Y  dx'j 


1/2 
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Figure  3:  Polynomial  Fitting  Of  The  Transformation  Curve 


(31) 


where  r(a:)  is  the  Gamma  function.  Apparently,  the  expected  magnitude  of  ^  will 
decrease  monotonically  as  long  as  the  window  widthcr  is  much  larger  than  one  pixel 
size,  which  is  usually  satisfied. 

Solving  Eq.  28  is  not  as  trivial  as  solving  Eq.  29.  But  since  the  expected  magnitude 
of  Uo  is  larger  than  those  of  Ui  and  U2,  we  propose  to  use  Eq.  28  as  an  iterative 
procedure  to  find  the  qth  time  estimation  of  such  as: 


(7, 


2(9) 


1  ln(f7o  +  fo(^o 

Jo 


^  ° - ^-U2)  -  In  Vo 


And  the  sufficient  condition  for  the  iteration  to  converge  is 

2  /oCl,  -  1(1  - 


r'(<T?>)  11=1 


fS  u„  + 


l<  1. 


(32) 


(33) 


For  some  locations,  the  convergence  condition  is  not  satisfied.  And  we  can  think  of 
the  convergence  condition  as  a  generalized  stability  criterion.  In  fact,  if  we  consider 
only  the  linear  term  in  the  Taylor’s  expansion,  this  convergence  condition  is  the  same 
as  the  stability  criterion  in  [6]. 
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3.2.3  Shift  Variance  Problem 


When  /i  7^  0  in  Eq.  22  and  Eq.  23,  we  must  consider  the  effect  of  the  shift  variance. 
Let  us  define, 


Vi  = 

r+oo 

J  —  OO 

ii[x)mi{x)dx , 

(34) 

r+oo 

/  p+oo  T  \ 

1  /  io(t)— 2^2  dt )  mi[x)dx. 

\J-oo  y2'Ks  J 

Vi  = 

J  —  OO 

(35) 

Erom  convolutions  of  the  two  images  with  moment  filters,  we  obtain  Ui  and  Vi. 
The  defined  Vi  can  not  be  obtained  from  convolutions.  But  to  estimate  s,  we  need  to 
compute  Vi  first.  In  other  words,  we  need  to  virtually  rotate  back  the  slanted  surface 
to  make  it  front-parallel,  then  compute  the  convolutions  Vi  of  the  rotated  surface. 
Replacing  Eq.  23  into  Eq.  34,  we  have 


mi{x)dx  -|- 


2^2 


xmi{x)dx. 


(36) 


Using  the  recursive  property  of  the  moment  filter  in  the  spatial  domain  in  Eq.  16, 
and  converting  the  above  equation  into  the  Fourier  domain,  we  obtain 

J  —  OO 

=  Vi  +  j5^^(cr^(Vi+3  +  2/oVj+2  +  /o  Vi+l)  — 

i(Vi+i-b2/oV  +  /jH_a)).  (37) 


In  other  words,  the  effect  of  the  linear  shift  variance  is  that  the  original  filter  outputs 
Vi  when  the  parameter  is  constant  within  the  window  are  recombined  linearly  to  form 

u. 

If  the  magnitude  of  ^  is  small,  we  can  then  ignore  the  terms  of  second  or  higher 
order  of  /x,  and  have  the  following  approximation  from  Eq.  37: 

V,  Vi-js^fi{a^{Vi+s  +  2foVi+2  +  f^Vi+^)- 
*(U+i  +  2/oU  +  /oK'-i) 

=  u  +  (38) 


i.e.  when  fi  is  small,  we  can  reverse  the  effect  of  the  linear  shift  variance,  and  obtain 
the  original  filter  output  Vi  when  the  the  parameter  is  a  constant  within  the  window. 
Once  we  have  computed  Vi,  we  can  compute  s  by  the  iteration  proposed  in  the 
previous  section  as 

=  ^  (\niUo  +  fos^U,  -  -  ln{Vo  +  ^V^o))  •  (39) 
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Now,  the  problem  is  that  we  don’t  know  fi.  Therefore,  we  also  need  an  algorithm 
to  estimate  Fortunately,  we  can  compute  ^  by  taking  the  ratio  of  the  zeroth  and 
the  first  moment  by  using  Eq.  27  and  Eq.  38: 

_  Ui  +  hs'‘Vi  -  p"(i  -  Us-‘)u^ 

Vo  +  aVa  ~  £/o  +  fosVi  -  p"(l  -  fSs^)U2  '  ' 

Combining  Eq.  39  and  Eq.  40  together,  we  have  two  equations  for  two  unknowns, 
i.e.  s  and  //.  Regarding  all  the  moments  as  constants  in  Eq.  39,  and  assuming  the 
iteration  in  Eq.  32  indeed  converges,  we  can  simplify  Eq.  39  as: 

s^  =  M{s^n),  (41) 

where  M.  can  be  regarded  as  a  complicated  function  of  /r.  Similarly,  Eq.  40  can  be 
simplified  as: 

(42) 

where  AT  is  a  complicated  function  of  s^. 

Eq.  41  and  Eq.  42  suggest  another  iterative  scheme  to  solve  both  s  and  ^  simulta¬ 
neously.  Let  u  =  and  v  =  and  suppose  the  9th  time  estimation  of  u  and  v  are 
and  respectively,  we  can  approach  the  true  value  of  u  and  v  by  iteration  as: 

(43) 

Z=  (44) 

And  the  sufficient  condition  for  the  iteration  to  converge  is  that  the  magnitudes  of 
both  the  complex  eigenvalues  Ei,i  =  0, 1  of  the  Jacobian  matrix  J  are  less  than  one, 

i.e. 

||Ei  ||<1.0,*  =  0,1,  (45) 

where, 

(46) 

Therefore,  when  the  convergence  condition  Eq.  45  is  satisfied,  we  can  iteratively 
estimate  depth  and  slope  simultaneously  as  in  Eq.  43  and  Eq.  44. 

3.2.4  Error  Estimation 

Error  estimation  is  necessary  if  we  have  to  merge  results  from  different  frequency 
bands.  The  convergence  conditions  (Eq.  33  and  Eq.  45)  can  identify  serious  cases  of 
the  weak  texture.  But  for  those  less  serious  cases,  we  still  need  to  quantify  the  depen¬ 
dency  of  errors  on  the  image  contents  in  order  to  optimally  merge  results  computed 
from  different  frequency  bands. 

Before  we  jump  into  the  estimation  of  errors,  we  need  to  be  clear  in  math.  Gen¬ 
erally,  both  Eq.  39  and  Eq.  40  will  return  complex  numbers  for  u  and  u,  which  in 
practice,  need  to  be  real  numbers.  Because  we  know  that  the  imaginary  part  of 


J  = 


0 

dM 

L 


dM 

0 
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computed  u  and  v  should  be  zero,  we  simply  drop  the  computed  imaginary  parts, 
and  therefore  keep  the  u  and  v  meaningful.  As  we  are  concerned  about  the  error 
estimation  of  u  and  n,  we  will  encounter  the  following  situation. 

Let  6  =  Sj.  +  j Si  as  a  zero-mean  complex  error  source,  whose  magnitude  E[SS*]  = 
m\  +  E\Sl\  =  2E161]  (we  assume  the  error  source  is  non-directional),  and 

u  =  AS 
V  —  BS^ 


we  then  have, 

E[nrVr\  =  El^~^  =  Re(.4i?*)i;[«>]/2.  (47) 

Generally,  if  u,  v  and  S  are  vectors,  and  A  and  B  are  matrices,  we  have, 

E[u,v^]  =  Re(A.E[.5.5*^]B*^)/2.  (48) 


There  are  two  types  of  error  sources  we  need  to  take  into  considerations.  The  first 
is  the  noises  in  images,  and  the  other  is  the  truncation  error  in  Eq.  26  and  Eq.  23. 
To  simplify  the  problem,  we  make  the  following  assumptions: 

1.  The  error  caused  by  truncation  in  Eq.  23  is  negligible. 

2.  The  noise  is  white  additive  noise. 


Based  on  these  two  types  of  error  sources,  the  error  in  estimating  u  and  v  in  Eq.  43 
and  Eq.  44  can  be  represented  as  (Appendix  C) 

=  AW,  (49) 

where  A  is  a  2x6  complex  matrix,  and  W  is  the  6x1  error  vector.  Applying  Eq.  48 
to  Eq.  49,  we  have  the  covariance  of  Ur  and  Vr  as 

E 

where  the  covariance  matrix  E^[WW*^]  is  a  6x6  diagonal  real  matrix. 


durduy  duydvy 
dUydVy  dVydVy 


=  Re(A£:[WW^]A*^)/2 


(50) 


du 

dv 


3.3  Moment  Filters  For  Stereo 

3.3.1  Problem  Definition 

The  method  of  depth  from  stereo  is  to  obtain  depth  information  using  a  quantitative 
measure  of  the  local  shift  between  two  images.  Assuming  zo(aj)  and  ii(x)  are  the  two 
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images  in  the  spatial  domain,  and  /o(/)  and  Ii{f)  are  in  the  Fourier  domain,  we  have 
the  following  relations. 


ii{x)  -  io{x  +  Do),  (51) 

h{f)  =  /o(/)e"°”.  (52) 

The  shift  between  two  images  is  characterized  by  the  disparity  Dq.  The  above 
relation  in  the  Fourier  domain  is  valid  only  when  Do  is  a  constant  within  the  window. 
If  Dq  varies,  instead,  in  the  neighborhood  of  window  center  xq  (denoted  as  e(a;o)),  i.e. 

Do  =  D  +  iJ,x,x  e  e{xo),  (53) 


we  then  have 


io(x  +  Do)  ^  ioix  +  D)  +  fix-^io(x  +  D),x  e  e(a;o), 

ax 


(64) 


where  we  truncated  terms  of  higher  order  of  fi  based  on  the  assumption  that  fi  is 
small. 

We  need  to  solve  the  following  two  problems: 

•  Finite  Window  Problem:  Since  Eq.  52  is  based  upon  an  infinite  window,  how 
can  we  compute  Dq  with  high  precision  when  the  window  is  finite? 

•  Shift  Variance  Problem:  Since  Eq.  52  is  based  upon  the  assumption  that  Dq  is  a 
constant  within  the  window,  how  can  we  compute  D  when  Dq  varies  within  the 
window  as  in  Eq.  53? 

3.3.2  Finite  Window  Problem 

In  this  section,  we  limit  our  discussion  to  the  situation  when  /u  =  0,  i.e.  the  disparity 
Do  is  indeed  constant  within  the  window.  Let  us  define  Ui  and  Vi  as  results  from 
convolving  the  two  images  with  the  I'th  moment  filter,  i.e. 


Ui  = 


Vi  = 


/  +  00 

io{x)mi{x)dx 

-OO 

/+00 

Io{fo)Mi{f)df 

-OO 

/  +  00 

ii{x)mi(x)dx 

-OO 

/-{-OO 

/„(/)e'4«"M.(/)rf/ 

-OO 


(55) 


(56) 


In  the  neighborhood  of  /o,  i.e.  e(/o),  the  transformation  curve  can  be  expanded 
according  to  Taylor’s  expansion 


JfDo 


(57) 
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where  higher  order  expansions  are  truncated. 

Replacing  Eq.  57  into  Eq.  56  and  using  the  recursive  property  of  the  moment  filter 
in  the  Eourier  domain,  we  have 

Eo  =  -  jDoU,  -  ^U2).  (58) 

Therefore,  we  have  the  equation  to  compute  Dq  as 


Do  =  j  [\n{Uo  -  jDoU^  -  -fU^)  -  In  Voj  (59) 

Again,  as  we  showed  in  the  focus  problem,  the  difference  between  Eq.  59  and  Eq.  52 
is  that  the  moment  filter  approach  uses  a  second  order  polynomial  to  approximate  the 
transformation  curve  in  the  neighborhood  of  /o,  while  the  traditional  approach 
uses  a  constant  to  approximate  the  curve.  As  we  can  see,  only  when  Ui  =  U2  = 

0,  is  Eq.  59  the  same  as  Eq.  52.  Therefore,  when  the  magnitude  of  Ui  or  U2  is  not 
small,  using  Eq.  52  to  compute  Dq  will  result  into  large  errors. 

By  analyzing  this  way,  it  is  not  a  surprise  to  learn  that  the  stability  constraint 
used  by  David  Fleet  is  actually  the  magnitude  ratio  of  Ui  and  Uq  ,  i.e. 

ai  laa  ||£/.  II 

where  (p  is  the  phase,  and  a  is  the  magnitude.  Eq.  60  also  justified  the  so-called 
phase  singularity,  which  is  location  where  Ui{fo)  is  so  large  compared  with  Uo{fo) 
that  Eq.  52  isn’t  even  approximately  valid.  We  can  see  that  such  situation  is  exactly 
what  Eq.  59  characterized. 

While  David  Fleet  figured  out  the  constraint  to  find  the  locations  of  singularity 
introduced  by  a  finite  window,  he  didn’t  find  a  specific  way  to  correct  the  error  caused 
by  a  finite  window.  As  the  same  as  we  did  in  the  problem  of  depth  from  defocus,  we 
can  use  Eq.  59  as  an  iterative  procedure  to  find  ^th  time  estimation  of  Dq  as 

=  j-  |^ln(17o  -  ^f/2)  -  InEoj 

=  (61) 
And  the  sufficient  condition  for  the  iteration  to  converge  is 

II  IK  1.  (62) 


3.3.3  Relations  to  the  Phase-Based  Method 

Some  researchers  [8,  23,  27]  recently  advocated  a  phase-based  approach  to  the  stereo 
problem.  This  approach  simply  divides  the  phase  difference  by  the  peak  frequency 
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or  the  instantaneous  frequency  to  infer  the  disparity.  Fleet [8]  showed  that  when 
the  stability  criterion  is  small,  more  precise  results  can  be  produced  by  using  the 
instantaneous  frequency  rather  than  the  peak  frequency. 

Suppose  we  truncate  the  Taylor’s  expansion  after  the  linear  term,  then  Eq.  58 
becomes 

=  (63) 

Rewriting  the  above  equation,  we  then  have 

jfoDo  =  In  Vo  -  ln(f/o  -  jDoUi).  (64) 

Because  the  stability  constraint  guarantees  that  ||  Ui  11<C||  Uq  ||,  the  second  term  of 
the  right  side  in  the  above  equation  can  be  approximated  as, 

ln(C/o  -  jDoUr)  ^InUo-  jDo^.  (65) 

b'O 

Replacing  Eq.  65  into  Eq.  64,  after  some  manipulations,  we  obtain, 

{fo-^)Do^j\nUo-j\nVo.  (66) 

Taking  the  real  part  of  the  both  sides  in  the  above  equation,  and  referring  to  the 
definition  of  the  instantaneous  frequency  in  Eq.  18,  we  can  see  that  the  above  equation 
literally  becomes  that  Do  multiplied  by  the  instantaneous  frequency  equals  the  phase 
difference.  Therefore,  we  can  regard  the  phase-based  approach  as  a  special  case  of 
the  moment  filter  approach  we  propose  here. 


3.3.4  Shift  Variance  Problem 


When  //  7^  0  in  Eq.  53  and  Eq.  54,  we  must  consider  the  effect  of  the  shift  variance. 
Let 


/■4-00 

Vi  =  ii{x)mi{x)dx , 

J  —  CO 

(67) 

/•+00 

El  =  /  io{x  -f  D)mi{x)dx 

J  —  OO 

(68) 

r+co 

J  — OO 

(69) 

From  convolutions  of  the  two  images  with  moment  filters,  we  obtain  Ui  and  Vi- 
The  define  Vi  can  not  be  computed  from  convolutions.  But  to  estimate  s,  we  need  to 
compute  Vi  first.  In  other  words,  we  need  to  virtually  rotate  back  the  slanted  surface 
to  make  it  front-parallel,  compute  the  convolution  of  the  rotated  surface.  Replacing 
Eq.  54  into  Eq.  67,  we  have 


/+00 

io{x  +  D)mi{x)dx  -f- 

-OO 

dio{x  +  D) 

/  - - xmAxjdx, 

J —OO  dx 


(70) 
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Using  the  recursive  property  of  the  moment  filter,  after  some  manipulations,  we 
can  rewrite  Eq.  70  as 


K-  =  Vi  -  p(i(/„Vi_.  -  Vi)  -  ct"(/cV,+,  -  Vi+2)). 


(71) 


If  jj,  is  small,  we  can  then  ignore  the  second  or  higher  order  terms  when  we  reverse 
the  above  equation,  i.e.  we  can  have  the  following  approximation: 


Vi  «  U  +  MW--1-^0-^'(/oK-+1-K+2)) 

=  U  +  flPi, 

where  Pi  is  used  for  simplicity  of  illustration. 

Combining  Eq.  59  and  Eq.  72,  we  have  the  equations  to  solve  D  and  yU, 


(72) 


D  = 


fo 


(ir{Uo  -  jDU^  -  -  ln(Uo  +  f^Po) 

U1+//P1  Ur-jDU2-lD^U^ 


(73) 


(74) 


Vo  +  fxPo  Uo-jDUi-^DW2 

The  same  iterative  scheme  we  used  in  focus  to  solve  s  and  yu  can  also  be  applied 
to  the  above  equations  to  solve  D  and  simultaneously. 

There  are  two  types  of  error  sources  that  we  need  to  take  into  considerations,  in 
error  estimation  as  we  did  in  focus.  The  first  one  is  noises  in  images,  and  the  other 
is  the  truncation  error  in  Taylor’s  expansion.  To  simplify  the  problem,  we  make  the 
following  assumptions  : 

1.  The  error  caused  by  truncation  in  Eq.  54  is  negligible. 

2.  The  noise  is  white  additive  noise. 


Based  on  these  two  types  of  error  sources,  the  error  in  estimating  D  and  fx  can  be 
represented  as  (Appendix  C) 

[  dD 

dfi 

where  A  is  a  2x6  complex  matrix,  and  W  is  a  6x1  error  vector.  Applying  Eq.  48  to 
Eq.  75,  we  have  the  covariance  of  Dr  and  fir  as 

^  dDrdDr  dDrdflr 
dDrdflr  dflrdflr 

where  the  covariance  matrix  E[WW*^]  is  a  6x6  diagonal  real  matrix. 


=  Real(AE[WW*^]A*^)/2,  (76) 


AW, 


(75) 
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3.4  Summary 

We  have  shown  how  the  moment  filters  can  be  used  solve  the  problems  of  focus  and 
stereo.  In  fact,  our  moment  filter  approach  can  handle  a  broad  category  of  problems, 
which  can  be  characterized  as  finding  a  single  parameter  in  the  transformation  be¬ 
tween  two  images.  As  long  as  the  transformation  is  well  modeled,  the  parameter  and 
its  gradient  can  be  accurately  estimated  from  one  frequency  band  by  the  following 
generic  algorithm: 

1.  In  the  neighborhood  of  the  peak  frequency  /o,  the  transformation  in  the  Fourier 
domain  can  always  be  approximated  as  a  polynomial  of  (/  —  /o),  using  Taylor’s 
expansion. 

2.  The  recnrsive  property  of  the  moment  filters  in  the  Fourier  domain  allows  ex¬ 
panding  the  convolution  of  the  Gabor  filter  with  one  image  to  be  expressed  as 
the  sum  of  convolutions  of  moment  filters  with  the  other  image.  Therefore,  we 
obtained  an  equation  of  the  parameter,  which  can  be  usnally  solved  by  iteration. 

3.  The  recursive  property  of  the  moment  filters  in  the  spatial  domain  allows  the 
modification  of  the  filter  outputs  caused  by  linear  shift  variance  of  the  parameter 
within  the  window  to  be  approximated  as  a  linear  combination  of  convolutions 
of  moment  filters.  This  modification  can  be  used  to  compute  the  gradient  of 
the  parameter  and  correct  the  error  in  estimating  the  parameter  caused  by  the 
gradient. 

Though  we  are  very  successful  in  solving  transformations  with  one  unknown  pa¬ 
rameter,  such  as  problems  of  focus  and  stereo,  there  are  some  intrinsic  limitations 
with  the  paradigm  based  on  one  frequency  band. 

•  It  is  usually  difficult  to  solve  transformations  which  have  more  than  one  param¬ 
eter  because  it  is  an  under-constrained  problem  for  one  frequency  band.  For 
example,  it  is  not  easy  to  extend  the  algorithm  of  stereo  to  solve  the  2D  cor¬ 
respondence  problem  because  2D  shift  involves  two  independent  parameters  in 
the  transformation.  Yet  because  we  can  have  tens  or  even  hundreds  of  frequency 
bands,  the  overall  problem  should  be  overconstrained  on  the  other  hand.  There¬ 
fore,  to  solve  problems  with  more  than  one  parameter,  we  can  no  longer  consider 
each  frequency  band  separately,  we  have  to  consider  all  together. 

•  There  is  another  related  problem,  which  we  refer  to  as  the  frequency  sampling 
problem.  Though  we  can  optimally  combine  results  from  different  frequency 
bands  through  Kalman  filtering,  the  selections  of  peak  frequency  fo  are  still 
ad  hoc.  They  by  no  means  represent  the  image  content  either  completely  or 
nonredundantly. 

To  overcome  these  problems,  we  propose  another  set  of  filters,  “hypergeometric” 
filters,  which  have  similar  recursive  properties  as  the  moment  filters  do,  and  don’t 
suffer  from  the  above  two  problems. 
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4  Hypergeometric  Filter  Approach 

Hypergeometric  filters  have  the  similar  recursive  properties  as  the  moment  filters 
have,  and  two  additional  properties: 

•  All  hypergeometric  filters  have  only  one  single  peak  frequency.  Because  the  high 
order  moment  filters  have  two  peak  frequencies,  we  have  to  base  our  analysis 
on  zero-order  filters  which  have  one  single  peak  frequency.  To  hypergeometric 
filters,  We  can  certainly  use  all  high  order  ones. 

•  The  set  of  hypergeometric  filters  samples  the  frequency  domain  completely  and 
nonredundantly. 

We  will  show  that,  by  using  the  hypergeometric  filters,  a  general  multiple  parameter 
extraction  problem  can  be  formulated  as  a  multidimensional  minimization  problem. 
Finally,  we  will  apply  this  new  technique  to  the  optical  flow  problem. 


4.1  Theory  of  Hypergeometric  Filters 

We  define  the  “hypergeometric  filters”  or  “H  filters”  in  the  Fourier  domain  as  (m  = 
l,2,...,oo): 


= 

Hoif)  = 


f  when  /  >  0 

I  0  when  /  <  0 

f  0  when  /  >  0 


where  Cm  is  a  real  normalization  constant  as 


=  2cr" 


TTcr 


r(m-M/2)' 


(77) 

(78) 

(79) 


(80) 


Intuitively,  a  H  filter  represents  either  one  peak  of  the  two  of  a  moment  filter  in 
Figure  2  with  fo  =  0.  Therefore,  we  can  easily  verify  that  Hm{f)  is  an  asymmetric 
bell-shaped  band  pass  filter  with  the  peak  frequency  at: 


Peak[i7,„(/)]  = 

(7 

(81) 

Peak[.ffo(/)]  = 

0, 

(82) 

Peak[i7_„,(/)]  = 

y/m 

(83) 

Doing  inverse  Fourier  transform  in  Hm{f)  in  Eq.  77  and  H-m{f)  in  Eq.  79,  we 
then  have  the  H  filters  represented  in  the  spatial  domain 


^m(^)  — 


7rr(m  -f  l/2)cr 


20-2 
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ZbI  order 


Sth  order 


lOUi  order 


Figure  4:  H  Filters  and  Their  Closest  Gabor  Filters 
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2  ’2’2<72b 


(84) 

(85) 
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where  $(ci,  c,  a;)  is  a  special  type  of  hypergeometric  function,  usually  called  “confluent 
hypergeometric  function”  or  “Rummer’s  function”  in  literature,  and  and  bm  are 
constants  decided  by  m.  Eq.  84  has  already  been  normalized  so  that  the  energy  of 
the  Alter  hm{x)  is  one.  The  solid  curves  in  Figure  4  are  H  filters  of  second,  fifth  and 
tenth  order. 

As  shown  in  Appendix  D,  the  H  filters  have  the  following  property: 


h2i{x)  +  h-2i{x) 


h2i-\-l{x)  —  h_^2i-\-l){x) 


jk2xL{i, 


1 


2’2(r2 


)e' 


20-2  , 


(87) 

(88) 


where  ki  and  k2  are  real  constants  decided  by  i,  and  L{n,a,x)  is  the  generalized 
Laguerre’s  polynomial.  Comparing  Eq.  87  and  Eq.  88  with  Eq.  2  and  Eq.  13,  referring 
to  the  fact  that  moment  filters  are  orthogonal  and  complete,  we  conclude  that  we  can 
losslessly  reconstruct  the  original  signal  from  coefficients  resulted  from  convolving  the 
image  with  H  filters. 

Another  important  property  of  H  filters  is  their  recursive  property  (Appendix  D) 
in  the  spatial  and  Fourier  domains,  i.e.  in  the  Fourier  domain, 

fHUf) 

fHoif) 


=  —Hm+lif), 

^m+1 

=  ^(7/,(/)  -  ff_i(/)). 

^m+1 
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(89) 

(90) 

(91) 


and  in  the  spatial  domain, 


xho{x) 

^  m  (^) 


-ja  ^^m  +  ^hra+ii^)-  ^ 
-^j(T{hi{x)  -  h_i{x)), 


m 


1/2 


1(^)  1  ? 


ja  I  \l m  +  ^h_(m+i)ix) 


m^ 


"'ij  m  —  1/2 


(m— 1)(^)  1  * 


(92) 

(93) 

(94) 


Unlike  Gabor  filters  which  have  positive  value  at  zero  frequency,  the  DC  compo¬ 
nents  of  H  filters  are  zero  except  for  m  =  0.  Therefore,  the  H  filters  have  the  advantage 
of  avoiding  the  strong  DC  bias  existing  in  Gabor  filters  which  are  extremely  harmful 
when  the  peak  frequency  is  low. 

Now  let’s  look  at  the  effective  bandwidth  in  the  Fourier  domain  and  effective  spatial 
extension.  According  to  Gabor  [9],  we  have 

n^x-‘  II  fe„(x)  IP  dx 
as  II  A„(x)  IP  dx  ’ 

-  Peak[g„(/)])^  II  g„(/)  IP  if 

SIS  II  HM)  IP  df 


^xl  = 

Ml  = 


Figure  5  shows  the  effective  bandwidth  and  effective  spatial  extension.  As  we  all 
know  the  uncertainty  principle  in  local  frequency  analysis,  the  effective  bandwidth 
and  effective  spatial  extension  must  satisfy  the  following  uncertainty  relation: 

=  sjAxlAfl  >\.  (97) 

And  Figure  6  shows  that,  as  Gabor  filters  have  the  minimal  uncertainty,  the  product 
of  the  spatial  and  the  spectral  uncertainties  of  the  H  filters  is  almost  the  minimal 
though  it  is  slightly  larger  (£„  <  0.51,  when  m  >  5). 

Since  Gabor  filters  are  the  only  filters  which  have  minimal  uncertainty,  and  the 
uncertainty  of  H  filters  approaches  the  minimal  value  as  m  increases,  we  can  predict 
that  the  H  filters  become  more  and  more  similar  to  Gabor  filters  as  m  increases.  The 
prediction  is  verified  by  Figure  4,  in  which  solid  curves  are  H  filters,  dashed  curves 
are  the  closest  Gabor  filters,  in  the  sense  that  they  have  the  same  peak  frequency  and 
the  same  bandwidth. 

The  H  filters  sampled  the  frequency  space  automatically  as  illustrated  in  Eq.  81 
and  Eq.  83.  Unlike  the  frequency  sampling  in  short  time  Fourier  transform,  which  is 
uniform,  the  H  filters  sample  the  frequency  space  in  a  non-uniform  fashion.  Figure  7 
shows  the  samplings  of  the  frequency  space  in  STFT  and  H  filters,  in  case  of  cr  =  2.0 
and  window  size  is  21.  Such  a  sampling  in  frequency  domain  will  be  shown  to  have 
great  advantage  in  combining  information  contained  in  different  frequency  bands  and 
solving  the  multiple  parameter  problem. 
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A  f  2 


Effective  Spatial  Extension 
Figure  5:  Uncertainty  in  The  Spatial  and  Frequency  Domains 


Figure  6:  Product  of  the  Spatial  and  Frequency  Uncertainty 


Figure  7:  Frequency  Sampling  in  STFT  and  H  Filters 


4.1.1  A  General  Approach  to  Multiple  Parameter  Problems 


As  an  abstract  example,  let  us  consider  a  problem  of  trying  to  recover  two  parameters 
ipiiPi)  ill  transformation  of  T(f‘,pi,p2)  between  two  images  Ii{f)  and  hif)- 
As  we  stated  before,  analyzing  a  single  frequency  band  as  we  did  in  the  moment 
filter  approach  is  not  enough  to  compute  two  parameters.  Therefore,  we  have  to  use 
multiple  frequency  bands  together  to  solve  this  problem.  For  simplicity  of  notation, 
we  only  use  H  filters  with  positive  frequency  peaks  hm{x),  (m  =  1, 2, ...)  in  this  section. 

The  transformation  can  be  represented  as 


W)  =  nf;puP2)W)  (98) 

Suppose  the  transformation  T{f-,pi,p2)  can  be  approximated  by  an  Nth  order  poly¬ 
nomial  Tp  in  the  neighborhood  of  Peak[i7„i(/)]  = 


N 


and  let 


Tpif-,Pi,P2)  =  1]  Ai(m,pi,p2)/' ,  /  e  e(^) 
i=o  ^ 


Kr, 


/  +  0O 

ii(x)hm(-x)dx 

-OO 

Jo 

/  +  00 

i2{x)hyn{-x)dx 

-OO 

=  h{f)Crare-^"’^"l^df 

Jo 

r+oo 

=  /  Hf)T{hp,,P2)<^r 

Jo 

and  use  the  recursive  property  of  the  H  filters  in  the  Fourier  domain. 


(99) 


(100) 


(101) 

(102) 
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we  then  have, 

^  c 

Vm=Yl  (103) 

j_0  Cm+i 

Apparently,  we  can  not  robustly  solve  (pi,p2)  from  a  single  equation  as  Eq.  103 
as  we  have  two  unknowns.  But  because  the  same  relation  should  be  satisfied  by  all 
frequency  bands,  we  can  therefore  minimize  the  following  to  find  (pi,P2): 

M2  N 

S{Pl,P2)=  II  Kn  -  I|^  (104) 

m=Mj  t=0 

where  Mi  and  M2  are  the  indexes  of  filters  which  have  the  lowest  and  highest  peak 
frequencies  respectively.  A  2D  conjugate  gradient  minimization  algorithm  can  be 
applied  to  the  above  function  to  detect  the  minima  quickly. 

We  could  also  represent  the  above  minimization  on  matrix  format.  Let, 

—^Ai{m,pi,p2)  =  kmi{pi,P2)-,  (105) 

we  can  rewrite  Eq.  103  as, 

Vi  '  ’  kiQ  kii  kif]  0  0  0  •  '  Ui' 

V2  0  k2Q  ^21  ■  ■  ■  ^2iV  0  0  •  •  •  U2 

Vz  =  0  0  ^30  ^31  •••  0  •••  Uz  '  (106) 

;  J  [o'-.  •••  •••  •••  •••  •••  0  J  L  . 

or  in  a  more  concise  form, 

V  =  K(;,„P2)U,  (107) 

where  V  and  U  are  known  vectors  from  convolutions,  K(pi,p2)  is  a  matrix  whose 
elements  are  functions  of  (pi,P2)-  Therefore,  the  minimization  of  S  in  Eq.  104  is 
equivalent  to  minimizing  the  norm  of  difference  between  the  left  and  the  right  side  of 
Eq.  107. 

Whenever  we  try  to  extract  multiple  parameters,  there  is  a  so-called  aperture 
problem,  which  is  caused  by  lack  of  constraint  along  certain  direction  in  the  parameter 
space.  In  context  of  multidimensional  minimization,  the  aperture  problem  happens 
when  the  function  S  has  a  canyon  instead  of  an  obvious  minima  in  the  parameter 
space.  Therefore,  the  aperture  problem  can  be  represented  by  an  error  covariance 
matrix  with  a  large  eigenvalue  along  one  direction  as  we  will  show  later  in  experiments. 

4.2  Hypergeometric  Filter  Approach  For  Optical  Flow 

We  will  apply  the  technique  of  solving  a  general  multiple  parameter  problem  to  the 
two-parameter  optical  flow  problem.  The  2D  transformation  in  optical  flow  problem  is 
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assumed  to  be  translation,  i.e.  assuming  two  adjacent  frames  are  io(x,y)  and 
and  the  image  velocity  is  we  then  have, 

ii{x,y)  =  ioix  +  v^,y  +  Vy),  (108) 

Hf.J.)  =  (109) 

For  simplicity  of  illustration,  we  only  approximate  the  above  transformation  by  a 
second-order  polynomial  according  to  Eq.  99  as 

rp  (  f  f  \  _  piifxV^+fyVy) 

^p\Jxi  Jyi  ^X1 '-'y)  — 

~  (1  JVxi^fx  fox)  "b  jf^yify  foy)  “  ^xifx  fox) 

-  MU,  -  M  - 

=  EEC’., (110) 

i=0  j=0 

Considering  only  positive  peak  frequency  H  filters,  i.e. 

Hmnifxjy)  =  HmUx)Hn{fy), 


we  then  obtain  the  following  relation  from  Eq.  103; 


T/  _  \  ^rn  Cn 

ymn  —  /  ,  /  ,  ^i.i\  )  I'^xi'^y)  ^(m+i)(n+j)' 

i=:0  j=0  ^  ^  Cn+j 

Therefore,  the  function  to  be  minimized  will  be. 


(111) 


^2  ^2  2  2-i  r  r 

SM  v,)=  ■£  j:  II  c«( v-  it  ■ 

m=Min=Ni  i=0  j=0  ^  ^  Cm+i^^n+j 

(112) 

H  filters  with  negative  or  zero  peak  frequency  can  be  similarly  applied,  and  they  can 
all  summed  into  S{vx,Vy). 


4.2.1  Error  Estimation 


Two  error  sources  are  going  to  be  considered  in  error  estimation.  One  is  the  white 
noise  contained  in  the  original  images.  The  other  is  the  truncation  error  in  Eq.  110. 
Suppose  the  intensity  of  the  white  noise  contained  in  the  two  images  are  both  it  is 
equivalent  to  that  the  first  image  contains  2N^  noise,  and  the  second  image  contains 
no  noise.  Therefore,  both  error  sources  result  in  inaccuracy  of  Vmn  in  Eq.  112,  which 
consequently  causes  the  estimation  error  of  (vx,Vy). 

Suppose  the  errors  of  Vmm  Vx  and  Vy  are  dVmn,  dvx  and  dvy  respectively,  we  can 
differentiate  both  sides  of  Eq.  112  as 


dvx 


M2  N2  2  2-i 

'Me  ■£  E((k:,-eec; 

m—Mi  n=Nl  i=Q  j=0 
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2  2-i  „  „ 

i=0  7=0  Cm+iC„+j 


t^(77l+2)(n+j)))7 


(113) 


Q  C  Ma  iVa  2  2—7 
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^(m+7)(77+j)) 


2  2-7  afi__ 

(  I'Tnt'Ti  Tr  X7 

_  „  ^(ni+7)(n+i)))- 
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(114) 


At  the  minima,  the  above  gradient  should  be  zero.  Then  the  task  of  the  error  es¬ 
timation  is  to  estimate  how  the  minima  will  move  when  there  is  a  small  random 
perturbation  in  Ktxti-  Let 


n  \  r<  i  \/^  Cm,Cn  jj 

^mnKPxi'^y)  /  .  /  .  i  I'^xi'^y)  ^ (tn-^i)(n+i)i 

7=0  j=0  ^  ^  Cm+iC„4.j 

we  can  take  differentiation  of  right  sides  of  Eq.  113  and  Eq.  114, 


(115) 
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Re  E(k:»  -c;.)- 
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dv,  j  ’ 


or  in  a  more  concise  form. 


a(  t‘ 


=  w, 
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where  A  is  a  2x2  real  matrix,  and  W  is  a  2x1  real  vector.  From  above  equation,  we 
can  have  the  covariance  error  estimation, 
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where  we  have 
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by  applying  Eq.  48. 
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4.3  Summary 

Our  new  hypergeometric  filter  approach  can  in  principle  handle  any  problem  of  ex¬ 
tracting  one  or  more  parameters  from  a  nonstationary  transform  between  two  images. 
We  showed  how  to  approach  the  optical  flow  computation  problem,  which  is  a  two 
parameter  extraction  problem.  For  different  parameter  extraction  problems,  the  al¬ 
gorithms  of  the  hypergeometric  approach  differ  only  in  the  Taylor’s  expansion  for 
different  transforms. 


5  Comparisons  of  Moment  Filters  and  Hyperge¬ 
ometric  Filters 

Because  both  moment  filters  and  hypergeometric  filters  have  similar  recursive  prop¬ 
erties  in  the  spatial  and  Fourier  domains,  they  have  the  same  strategy  for  solving 
parameter  extraction  problems.  This  strategy  uses  a  polynomial  to  fit  the  underlying 
transformation  curve,  which  is  a  surface  in  two  parameter  cases  and  hypersurface  in 
more  than  two  parameter  cases,  and  therefore  obtain  results  of  much  higher  precision 
the  traditional  approach  which  uses  a  constant  to  approximate  the  curve.  In  this  re¬ 
spect  the  two  approaches  are  equally  effective  in  computing  accurately  and  handling 
the  foreshortening  problem.  On  the  other  hand,  the  hypergeometric  filter  approach 
can  solve  multiple  parameter  problems,  while  the  moment  filter  is  good  at  dealing 
with  one  parameter  problems  even  though  it  is  possible  that  multiple  bands  can  be 
combined  to  solve  multiple  parameter  problems. 

The  moment  filter  approach  is  based  oh  a  single  frequency  band.  Thus  we  have  the 
freedom  to  sample  the  frequency  space  as  we  like.  The  disadvantages  of  this  freedom 
is  that  there  is  no  way  that  we  can  sample  the  frequency  domain  perfectly.  A  dense 
sampling  scheme  undermines  the  assumption  of  the  Kalman  filtering  in  combining  re¬ 
sults,  while  a  sparse  sampling  scheme  abandons  information  contained  in  images.  On 
the  other  hand,  the  advantage  of  this  freedom  is  that  we  can  have  more  sophisticated 
window  scheme  such  as  the  wavelet  scheme  in  [30]. 

The  advantage  and  disadvantage  of  the  hypergeometric  filter  approach  are  exactly 
the  opposite  of  those  of  the  moment  filter  approach.  The  hypergeometric  filter  ap¬ 
proach  samples  the  frequency  space  in  a  canonical  fashion  such  that  all  information 
are  utilized.  But  on  the  other  hand,  we  don’t  have  the  freedom  to  choose  different 
window  size  for  different  frequency  bands. 


6  Implementation  Issues  and  Experiments 

6.1  Moment  Filter  Approach 

In  this  section,  we  will  explain  some  implementational  issues  of  the  moment  filter  ap¬ 
proach,  and  quantitatively  compare  the  focus  and  stereo  algorithms  based  on  moment 
filters  with  other  algorithms  in  literature. 
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Figure  8;  Frequency  Sampling  in  Fixed  Window  Scheme 

6.1.1  Sampling  of  the  Frequency  Domain 

The  moment  filter  approach  described  above  is  based  on  one  single  Gabor  filter.  Since, 
the  information  contained  in  a  image  is  usually  distributed  in  all  frequency  bands,  we 
have  to  use  multiple  Gabor  filters  and  their  accompanying  moment  filter  sets.  At  this 
point,  the  sampling  of  the  frequency  domain  by  Gabor  filters  is  still  arbitrary.  In  our 
2D  implementation  of  the  moment  filter  approach,  the  sampling  is  done  in  a  polar 
coordinate,  i.e.  the  frequency  domain  is  sampled  in  radial  and  angular  directions.  The 
angular  sampling  rate  is  one  per  30  degrees,  and  the  radial  sampling  rate  is  chosen 
as  one  per  0.7/<t,  where  cr  is  the  window  size  parameter.  The  minimum  frequency  is 
one-tenth  of  the  Nyquist  frequency,  and  we  use  only  the  first  48  samples  in  the  polar 
sampling.  Therefore,  we  have  totally  48  Gabor  filters  and  48  moment  filter  sets.  In 
the  experiments  shown  below,  we  use  up  to  the  4th  order  moment  filters. 

Though  the  moment  filter  approach  can  be  applied  under  any  frequency  sampling 
scheme,  we  have  to  choose  one  for  experimentation.  There  are  two  window  schemes 
which  are  related  to  the  radial  frequency  sampling.  In  the  fixed  window  scheme, 
because  cr  are  the  same  for  all  Gabor  filters,  the  radial  frequency  sampling  is  uniform. 
In  the  variable  window  scheme, 

<7  =  fcA,  (121) 

where  A:  is  a  constant  and  A  is  the  peak  wavelength  of  the  Gabor  filter,  the  radial 
frequency  sampling  is  nonuniform.  We  don’t  intend  to  compare  these  two  window 
schemes  here,  but  rather  provide  two  choices  for  the  convenience  of  experimentation. 
Figure  8  and  Figure  9  illustrate  the  sampling  of  the  frequency  plane  under  the  window 
schemes.  The  circles  represent  the  effective  bandwidth  of  Gabor  filters. 

The  major  difficulty  of  the  filter-based  approach  in  implementation  is  its  enor¬ 
mous  amount  of  computation  and  memory  requirements  if  we  need  to  sample  many 
frequency  bands.  We  believe  that  the  as  parallel  computers  become  more  powerful, 
this  problem  will  be  much  less  severe.  Because  of  this  limitation,  we  only  use  48 
filters  which  is  a  very  sparse  set  of  samples  in  the  frequency  domain.  We  believe  that 
if  we  sample  the  frequency  domain  in  a  denser  fashion,  we  could  be  able  to  further 
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Figure  9:  Frequency  Sampling  in  Variable  Window  Scheme 
improve  the  results  reported  here. 

6.1.2  Focus 

We  implement  the  moment  filter  approach  for  the  focus  problem  using  the  third  order 
Taylor’s  expansion  in  Eq.  26.  Because  computation  of  Pi  in  Eq.  40  requires  the  fourth 
order  moment  filter,  we  need  up  to  the  fourth  order  moment  filters  for  every  Gabor 
filter.  To  accommodate  numerical  errors,  the  convergence  thresholds  in  Eq.  33  and 
Eq.  45  are  set  to  0.9.  Results  from  different  frequency  bands  are  merged  optimally 
through  Kalman  filtering  as  in  [24]. 

Subbarao’s  algorithm  in  [26]  is  directly  from  Eq.  21,  which  can  be  regarded  as  a 
special  case  of  our  moment  filter  approach,  in  which  the  Taylor’s  expansion  in  Eq.  26 
only  keeps  the  constant  term.  But  Subbarao’s  algorithm  can  not  be  extended  to  deal 
with  the  shift  variance  problem. 

We  implemented  both  algorithms  under  the  variable  window  scheme,  in  which  k 
in  Eq.  121  is  0.8.  The  performance  tests  are  done  on  both  synthetic  and  real  images. 
The  original  image  is  shown  in  Figure  10.  We  uniformly  blur  the  image  by  a  Gaussian 
function  with  parameter  <To  =  1.0  to  get  the  first  synthetic  image  pair.  We  also  non- 
uniformly  blur  the  original  image  by  a  Gaussian  function  with  parameter  <to  increases 
linearly  as  the  column  increases,  to  get  the  second  synthetic  image  pair. 

Since  in  2D  images,  the  slope  of  a  surface  has  to  be  represented  by  a  vector  instead 
of  one  number  in  Eq.  22,  we  have  to  solve  three  unknown  equations  in  case  of  2D 
image  pairs  instead  of  two  unknown  equations  in  ID  case. 

We  compare  three  different  algorithms,  namely,  Subbarao’s  algorithm  (SA),  the 
moment  filter  algorithm  without  considering  slope  (MFFl),  and  the  moment  filter 
algorithm  with  considering  slope  (MFF2).  Table  1  shows  the  root  mean  square  errors 
of  the  synthetic  image  pairs  using  different  algorithms.  In  the  case  of  the  first  pair, 
MFFl  and  MFF2  shows  almost  no  error  at  all,  while  SA  has  a  significant  amount  of 
error.  In  the  case  of  the  second  pair,  the  RMS  error  of  SA  is  4  times  larger  than  that 
of  MFFl,  and  27  times  larger  than  that  of  MFF2!  We  can  see  that  the  main  error 
source  of  Subbarao’s  algorithm  is  from  the  truncation  of  information  from  moment 
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Uniformly  Blurred  Image  Nonuniformally  Blurred  Image 


Figure  10:  Synthetic  Texture  Images  For  Focus 


SA 

MFFl 

MFF2 

The  First  Pair 

0.070518 

0.000749 

0.000316 

The  Second  Pair 

0.190362 

0.046419 

0.007976 

Table  1:  RMS  Errors  of  Different  Algorithms 


filters.  When  the  surface  is  slanted,  i.e.  the  foreshortening  happens,  it  also  causes 
errors  depends  on  how  slanted  the  surface  is.  The  MFF2  algorithm  can  correct  the 
errors  caused  by  the  foreshortening. 

To  see  the  precision  of  different  algorithms.  Figure  11  shows  the  computed  local 
blurring  difference  ctq  of  the  first  image  pair,  and  Figure  12  shows  the  computed  local 
blurring  difference  <Tq  of  the  second  pair.  Furthermore,  Figure  13  illustrates  the  slope 
vectors  computed  by  MFF2  for  the  second  image  pair.  Note  that  the  slope  vector 
has  been  scaled  by  as  the  left  side  of  Eq.  42,  which  causes  the  magnitudes  of  slope 
vectors  increases  horizontally.  We  can  see  that  MFF2  not  only  can  correct  the  errors 
caused  by  the  foreshortening,  but  also  estimate  the  slope  of  the  surface. 

The  real  image  pair  is  shown  in  Figure  14.  These  two  pictures  are  taken  by  the 


Figure  11:  Computed  Blurring  Difference  of  The  First  Pair 
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SA 

Figure  12;  Computed  Blurring  Difference  of  The  Second  Pair 
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Figure  14:  Real  Image  Pair  of  A  Toy  House 


Photometries  camera  [28]  under  different  aperture  and  exposure  time  setting.  The 
toy  house  is  about  two  meters  away  from  the  camera,  and  the  size  of  the  house  itself 
is  about  2.5  inches  wide,  3.0  inches  long,  and  4.5  inches  high.  To  compensate  for 
the  optical  vignetting  and  different  exposure  amount  between  two  images,  we  use  an 
approximately  uniform  illumination  to  correct  them  by  putting  a  diffuser  in  front  of 
the  camera,  and  a  light  source  far  away. 

Figure  15  shows  the  computed  blurring  difference  of  the  house  image  pair  using 
Subbarao’s  algorithm.  We  can  see  that  the  result  is  so  noisy  that  the  shape  of  the 
house  can  not  be  observed  with  confidence.  Figure  16  shows  the  result  using  MFFl. 
By  now  the  shape  of  the  house  is  clearly  visible  though  the  discontinuity  still  generates 
a  lot  of  spikes.  Finally  Figure  17  shows  the  result  using  MFF2.  We  can  see  that  the 
result  is  even  more  smooth.  Additionally,  Figure  18  shows  the  computed  slope  and 
estimated  error  {E[du.rdur\  in  Eq.  50)  in  MFF2.  We  can  see  that  most  of  the  slope 
vectors  have  the  correct  direction  despite  that  the  slope  is  very  small  [pt  0.002  in 
Eq.  22),  and  the  discontinuity  generates  very  large  slope  vectors  and  estimated  errors 
which  can  be  used  to  locate  the  discontinuity  in  later  processing  stages. 

6.1.3  Stereo 

Though  filter-based  stereo  algorithms  may  be  able  to  alleviate  the  mismatch  problem, 
our  concern  here  is  focused  on  the  precision  of  subpixel  registration.  We  will  compare 
our  algorithms  with  Kanade  &  Okutomi’s  adaptive  window  scheme  [13].  A  simple 
SSD-based  matching  is  performed  to  generate  pixel  resolution  disparity,  and  the  same 
pixel  resolution  disparity  is  used  for  our  algorithms  and  Kanade  &  Okutomi’s  algo¬ 
rithm. 

We  implemented  out  algorithms  under  the  fixed  window  scheme,  in  which  the 
window  size  a  =  7.0.  The  minimal  radial  frequency  is  tt/IO,  and  the  maximal  radial 
frequency  is  about  tt/S.  There  are  totally  48  filters  as  shown  before.  For  Kanade  & 
Okutomi’s  algorithm,  the  range  of  the  window  size  is  from  3  to  21,  and  the  iteration 
is  5  times.  For  simplicity,  we  will  denote  Kanade  &  Okutomi’s  algorithm  as  KOA, 
the  moment  filter  approach  without  slope  correction  as  MFSl,  and  the  moment  filter 
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Figure  18:  Slope  and  Error  Estimations  From  MFF2 
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Figure  19;  The  Original  Texture  Image 

approach  with  slope  correction  as  MFS2.  The  RMS  errors  we  will  use  below  are  in 
pixel  width. 

We  first  test  stereo  algorithms  on  synthetic  images.  The  first  two  test  image  pairs 
are  obtained  by  artificially  shifting  an  original  image.  The  first  image  is  by  a  uniform 
subpixel  shift  of  the  original  one,  and  the  second  image  by  a  nonuniform  subpixel 
shift,  which  naturally  results  in  stretching  or  compression.  Technically  the  subpixel 
shift  is  achieved  by  linearly  interpolating  the  original  image.  Figure  19,  Figure  20  and 
Figure  21  shows  the  original  and  shifted  images.  The  uniform  shift  in  Figure  20  is  1.55 
pixel  size,  and  the  nonuniform  shift  in  Figure  21  increase  linearly  from  1.0  pixel  size 
at  zero  column  to  3.54  pixel  size  at  the  rightmost  column.  The  third  synthetic  image 
pair  is  the  20th  and  21st  frames  of  the  translating  tree  sequence  [3]  as  in  Figure  22. 
The  scene  is  a  planar  surface  textured  by  a  picture  of  a  tree.  Because  the  ground 
truth  of  the  motion  is  in  horizontal  direction,  we  can  use  them  as  a  stereo  pair.  The 
disparity  is  between  1.73  and  2.26,  and  it  increases  linearly  from  left  to  right.  The 
interpolation  method  for  subpixel  shifting  is  unknown. 

Table  2  shows  the  RMS  errors  of  disparity  values  computed  by  different  algorithms. 
We  see  that  for  the  first  uniformly  shifted  image  pair,  all  algorithms  can  do  very  well. 
Because  the  RMS  errors  in  this  case  are  less  than  one  percent  of  the  pixel  size, 
we  believe  that  we  are  approaching  the  limit  set  by  the  spatial  sampling.  For  the 
nonuniformly  shifted  pair,  we  see  that  the  errors  are  significantly  larger  than  those 
from  the  uniformly  shifted  pair.  Figure  23  shows  the  disparity  values  of  an  arbitrary 
row  for  the  nonuniformly  shifted  image  pair.  We  see  that  the  ripples  generated  by 
KOA  are  significant  in  this  example.  Furthermore,  Figure  24  illustrates  the  slope 
computed  by  MFS2.  For  the  tree  pair.  Figure  25  shows  the  disparity  maps  computed 
by  different  algorithms.  And  Figure  26  illustrates  the  slope  computed  by  MFS2. 

We  tested  these  algorithms  on  two  real  image  pairs,  which  were  taken  by  the 
Photometries  Camera  in  CIL[28].  The  targets  were  about  2.5  meters  from  the  camera, 
and  the  camera  was  translated  horizontally  by  0.5  inch  between  left  and  right  view. 


36 


Figure  20:  The  Uniformly  Shifted  Image 


Figure  21:  The  Nonuniformly  Shifted  Image 


Uniformly  Shifted  Pair 

Nonuniformly  Shifted  Pair 

Tree  Pair 

KOA 

0.008 

0.035 

0.016 

MFSl 

0.008 

0.024 

0.012 

MFS2 

0.006 

0.020 

0.010 

Table  2:  RMS  Errors  of  Disparity  From  Synthetic  Images 
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Figure  22:  The  20th  and  21st  Frames  From  Translating  Tree  Sequence 


Disparity 


Figure  23:  Computed  Disparity  From  the  Nonuniformly  Shifted  Pair 
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Top  Patch 
(6690  Pixels) 

Left  Patch 
(16025  Pixels) 

Right  Patch 
(19625  Pixels) 

KOA 

0.063 

0.032 

0.027 

MFSl 

0.039 

0.024 

0.020 

MFS2 

0.018 

0.020 

0.016 

Table  3;  RMS  Errors  of  The  Cube 


Figure  27:  Images  of  A  Textured  Cube 

Because  there  is  no  vergence,  the  disparity  map  of  a  planar  surface  should  still  be 
planar.  We  will  fit  a  plane  to  regions  in  the  disparity  map  which  correspond  to  planar 
surfaces,  and  then  compare  quantitatively  the  performance  of  stereo  algorithms.  Note 
that  when  we  fit  planar  patches,  we  exclude  those  pixel  near  boundary  so  that  the 
RMS  errors  we  computed  are  not  influenced  by  the  depth  discontinuities. 

The  first  pair  is  images  of  a  textured  cube  as  in  Figure  27.  There  are  three  large 
planar  surfaces,  which  we  denote  as  top  patch,  left  patch,  and  right  patch  as  in 
Figure  28.  The  difference  between  the  largest  and  the  smallest  disparity  is  about  3.0 
pixel  widths.  Figure  29,  Figure  30,  and  Figure  31  show  the  computed  disparity  maps. 
We  can  see  that  the  disparity  map  computed  from  KOA  contains  more  ripples  than 
those  from  MFSl  and  MFS2.  Furthermore,  Figure  32  illustrates  the  computed  slope 
from  MFS2,  which  correctly  indicates  not  only  the  direction  of  the  slope,  but  also  the 
magnitude  of  the  slantness.  Table  3  shows  the  RMS  error  after  plane  fittings.  We 
observe  that  both  MFSl  and  MFS  have  better  performance  than  KOA,  and  the  RMS 
error  of  MFS2  keeps  almost  constant  when  the  surface  slantness  increases,  which  we 
attribute  to  the  slope  correction  MFS2  performed. 

The  second  real  pair  is  images  of  the  same  toy  house  we  used  in  the  focus  experi¬ 
ments.  Figure  33  shows  the  image  pair  we  used  for  the  stereo  computation.  There  are 
two  slightly  slanted  planar  surfaces,  which  we  denote  as  the  left  patch  and  the  right 
patch  as  in  Figure  34.  The  difference  between  the  largest  and  the  smallest  disparity 
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Figure  29:  Disparity  Map  of  The  Cube  From  KOA 
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Figure  31:  Disparity  Map  of  The  Cube  From  MFS2 
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Figure  32;  Slope  of  The  Cube  From  MFS2 


Left  Patch 
(11260  Pixels) 

Right  Patch 
(13204  Pixels) 

KOA 

0.032 

0.031 

MFSl 

0.028 

0.017 

MFS2 

0.019 

0.017 

Table  4:  RMS  Errors  of  The  Toy  House 


is  only  about  0.5  pixel  size,  excluding  the  background.  Figure  35,  Figure  36  and 
Figure  37  show  the  computed  disparity  maps  from  different  algorithms.  Again,  we 
can  see  the  consistent  step-like  ripples  from  KOA.  Figure  38  illustrates  the  computed 
slope  from  MFS2,  which  indicates  the  and  magnitude  and  direction  of  the  slope  on 
the  surfaces,  and  locates  the  discontinuities  by  large  slope  values.  Table  4  shows  the 
RMS  error  after  plane  fittings  of  the  two  patches. 


6.2  Hypergeometric  Filter  Approach 

6.2.1  Computation  of  Hypergeometric  Filters 

Unlike  most  other  filters,  the  numerical  computation  of  a  hypergeometric  filter  itself 
is  a  nontrivial  problem.  The  infinite  series  representation  of  the  confluent  hypergeo¬ 
metric  function  in  [25]  (Equation  47:3:1)  converges  very  slowly  when  the  parameters 
are  not  very  small.  Furthermore,  when  one  parameter  is  negative,  the  adjacent  terms 
in  the  infinite  series  could  conceal  each  other  such  that  their  sum  is  dominated  by 
roundoff  errors.  From  experiments,  we  found  that  the  so-called  “Algorithm  707”  in 
[19,  18]  worked  well  most  of  time  in  computation  of  the  hypergeometric  filters  with 
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Figure  33:  Stereo  Images  of  A  Toy  House 


Figure  34:  Planar  Surfaces  of  The  Toy  House 
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Figure  38:  Slope  of  the  House  From  MFS2 
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IP  =  10.  But  it  did  fail  sometimes  indicated  by  either  the  numerical  evaluator  of  the 
confluent  hypergeometric  function  failed  to  terminate  or  the  computed  hypergeomet¬ 
ric  Alter  had  a  sharp  spike.  In  case  of  its  failures,  we  observed  that  applying  Kummer’s 
transformation  first,  and  then  evaluating  the  confluent  hypergeometric  function  can 
solve  this  difficulty.  Kummer’s  transformation  is 

$(a,c,  x)  =  e“’$(c  —  a,c,  — a;).  (122) 

Fortunately  we  don’t  have  to  generate  the  H  filters  online,  which  may  be  a  very 
challenging  numerical  problem.  Instead,  we  can  generate  the  H  filters  according 
to  Eq.  84,  and  in  case  the  numerical  algorithm  fails  for  a  filter,  we  use  Kummer’s 
transformation,  and  recompute  the  filter.  Our  experience  is  that,  when  m  <  40, 
we  should  apply  the  Kummer’s  transformation  first  to  Eq.  84,  otherwise,  we  should 
evaluate  Eq.  84  directly. 

The  2D  extension  of  the  H  filter  is  as  straightforward  as: 

hij{x,y)  =  hi{x)hj{y).  (123) 

Eor  real  value  images,  because  the  Fourier  transform  I{fx,fy)  has  the  property, 

lif.Jy)  =  -fy), 

only  half  of  the  plane  needs  to  be  sampled.  From  Eq.  81,  the  Nyquist  frequency  tt  is 
reached  when 

m  cr^TT.  (124) 

Due  to  hardware  limitations,  we  will  not  be  able  to  use  all  the  filters  up  to  the  Nyquist 
frequency  in  our  implementation  here.  Instead,  we  will  use  the  {i,j)  (Eq.  123), 
correspondingly  {fxify)i  contained  in  the  shadow  area  in  Eigure  39,  in  which. 


Ui  = 

VMi 

a 

(125) 

(126) 

U2  = 

a 

Vi  = 

a 

(127) 

V2  = 

y/N~2 

a 

(128) 

Again,  we  believe  that  once  we  are  able  to  fully  explore  the  whole  Eourier  domain 
up  to  the  Nyquist  frequency  by  advancement  of  more  powerful  hardware,  the  results 
reported  here  should  be  improved  considerably. 

6.2.2  Optical  Flow 

We  adopt  the  2D  conjugate  gradient  algorithm  in  [22]  (function  frprmn)  to  com¬ 
pute  flow  vectors.  Eor  stability  of  numerical  computation,  we  actually  minimize  the 
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Figure  39:  Frequency  Sampling  in  Our  Implementation 


normalized  version  of  S{vj;,Vy)  in  Eq.  112, 


S  {Vx ,  ^2/) 


S{Vx,Vy) 


Em  E^ 


V  2' 

^  mn 


(129) 


As  we  did  in  the  moment  filter  approach,  we  truncate  the  Taylor’s  expansion  after 
the  third  order  term  in  Eq.  110. 

There  have  been  many  algorithms  to  compute  optical  flow  from  two  or  more  images. 
Among  the  published  work,  Barron  et  al  [3]  provided  a  quantitative  measure  for 
a  few  of  the  optical  flow  techniques.  We  will  demonstrate  the  capability  of  our 
new  technique  on  the  same  test  images  for  comparison.  While  most  of  techniques 
presented  in  [3]  use  a  sequence  of  images,  our  technique  use  just  two  adjacent  frames. 
No  pre-smoothing  in  the  spatial  domain  or  the  temporal  domain  is  done,  therefore 
the  quantization  error  and  noise  is  much  higher.  The  error  measurement  we  use  is 
the  angular  measure  in  [3], 


=  arccos 


‘^x'^Ox  T  ^3/^Oy  T  I 


Ox  +  +  1 


(130) 


where  {vx,  Vy)  is  the  computed  optical  flow,  and  (uoai,  ?^oy)  is  the  true  optical  flow.  We 
don’t  use  the  standard  deviation  measure  in  [3]  because  it  is  hard  to  interpret.  All 
angular  errors  are  in  degree  (°). 

Without  any  post  processing,  the  approach  proposed  here  will  certainly  compute 
a  velocity  at  every  pixel.  Two  well-known  intrinsic  problems  have  to  be  addressed, 
i.e.  the  aperture  problem  and  the  no-texture  problem,  both  resulting  from  a  lack  of 
intensity  variation  around  a  pixel.  If  the  error  estimation  procedure  indeed  works 
properly,  the  covariance  matrix  it  computed  will  have  one  large  eigenvalue  in  case 
of  the  aperture  problem,  and  two  large  eigenvalues  in  case  of  no  texture  at  all.  The 
corresponding  eigenvector  should  indicate  the  direction  in  which  large  uncertainty 
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Technique 

Number 

Horn  &  Schunck  (original) 

1 

Horn  &  Schunck  (original)  ||  V/  >  5.0 

2 

Horn  &  Schunck  (modified) 

3 

Horn  &  Schunck  (modified)  ||  V/  >  5.0 

4 

Lucas  &  Kanade  (A2  >  1.0) 

5 

Lucas  &  Kanade  (A2  >  5.0) 

6 

Uras  et  al.  (unthresholded) 

7 

Uras  et  al.  (det(Lf)  >  1.0) 

8 

Nagel 

9 

Nagel  V7  2>  5.0 

10 

Anandan 

11 

Singh  (step  1,  n  =  2,rw  =  2) 

12 

Singh  (step  1,  n  =  2,ro  =  2,  Ai  <  5.0) 

13 

Singh  (step  2,  n  =  2,  ru  =  2) 

14 

Singh  (step  2,  n  =  2,  tn  =  2,  Ai  <  0.1) 

15 

Heeger 

16 

Heeger  (level  0) 

17 

Heeger  (level  1) 

18 

Heeger  (level  2) 

19 

Waxman  et  al.  (<t/  =  2.0) 

20 

Fleet  Sz  Jepson  (r  =  1.0) 

21 

Fleet  &  Jepson  (r  =  1.25) 

22 

Fleet  &  Jepson  (r  =  2.5) 

23 

Table  5:  Numbering  the  Optical  Flow  Techniques 


occurs.  To  compare  our  results  with  those  from  other  techniques,  we  threshold  the 
computed  optical  flow  by  the  following  formulae: 


Ema,x 
vl  +  vl  +  1 


<Te 


(131) 


in  which  Emax  is  the  largest  eigenvalue  of  the  covariance  matrix,  and  Te  is  the 
threshold.  By  setting  Te  to  different  values,  we  can  obtain  a  curve  of  error  versus 
density  as  we  will  show  later. 

For  the  convenience  of  comparisons,  we  number  the  various  techniques  imple¬ 
mented  in  [3]  as  in  Table  5. 

For  the  Yosemite  sequence  (Figure  40),  we  use  the  9th  and  10th  frames  to  compute 
the  optical  flow  as  they  are  the  two  frames  in  the  middle  of  the  sequence.  The  size 
of  the  H  Alters  is  ct  =  7.0,  and  the  sampling  limits  in  Figure  39  are:  Mi  =  Ni  = 
0,M2  =  N2  =  10,  which  correspond  to  the  area  in  the  Fourier  domain:  (0.0  < 
fx  <  0.45,-0.45  <  fy  <  0.45).  In  other  words,  we  only  use  information  whose 
frequency  is  lower  than  one-sixth  of  the  Nyquist  frequency.  Figure  41  shows  the  raw 
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The  Ninth  Frame  in  Yosemite  Sequence  The  Correct  Flow 


Figure  40:  Yosemite  Image  Sequence 

(unthresholded)  optical  flow  computed.  Surprisingly,  the  movement  of  the  clouds 
is  captured  accurately.  The  curve  of  the  density  versus  average  error  is  shown  in 
Figure  42  along  with  the  results  reported  in  [3].  For  further  comparison,  we  listed  the 
average  error  and  the  corresponding  density  in  Table  6.  The  numbers  in  Figure  42 
are  the  numbered  techniques  in  Table  5.  We  can  see  that  except  for  a  couple  of 
techniques  which  impose  smoothness  constraint  at  100%  density,  our  approach  has 
the  best  results  for  all  densities,  despite  the  fact  that  we  only  use  low  frequency 
information  from  only  two  frames. 

The  monotonic  increase  of  the  average  error  with  respect  to  more  restrictive  eigen¬ 
value  threshold  indicates  that  the  error  estimation  indeed  represented  the  uncertainty 
properly.  Furthermore,  as  a  covariance  matrix  can  be  represented  by  a  ellipse  whose 
two  principal  directions  are  the  two  eigenvectors,  and  whose  radii  are  the  square  roots 
of  the  corresponding  eigenvalues,  we  superimpose  the  ellipses  on  top  of  the  darkened 
image  in  Figure  43.  First  of  all,  from  the  size  of  ellipses,  it  indicates  that  the  top 
portion  and  low-left  portion  have  large  errors  which  is  consistent  with  Figure  41  and 
Figure  40.  Also  we  observe  that  at  locations  where  the  aperture  problem  is  obvious, 
e.g.  near  a  single  edge,  the  ellipses  are  elongated  in  the  direction  of  edges,  which 
represent  the  lack  of  information  in  those  directions. 

The  second  synthetic  sequence  is  the  translating-tree  sequence  (Figure  44)  from  [3]. 
We  used  the  20th  and  21st  frames  to  compute  the  optical  flow.  The  size  of  the  H  filters 
is  (T  =  3.5,  and  the  sampling  limits  in  Figure  39  are:  M\  =  N\  =  OjMa  =  N2  =  15, 
which  correspond  to  the  area  in  the  Fourier  domain  (0.0  <  /r  <  1.233,-1.233  < 
fy  <  1.233).  In  other  words,  the  highest  frequency  of  the  information  we  used  is 
one-third  of  the  Nyquist  frequency.  Figure  45  shows  the  unthresholded  optical  flow. 
Except  sparse  locations  where  there  are  no  texture,  the  computed  optical  flow  is 
pretty  consistent  with  the  ground  truth.  The  curve  of  the  density  versus  average 
error  is  shown  in  Figure  46,  in  which  some  results  with  large  errors  reported  in  [3] 
are  not  included,  and  and  Table  7.  Again  we  can  see  that  our  approach  based  upon 
only  two  adjacent  frames  and  only  low  frequency  information  outperformed  all  the 
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Average  Angular  Error 


Figure  41:  Computed  Optical  Flow  of  the  Yosemite  Scene 


density 

Figure  42:  Average  Error  For  the  Yosemite  Scene 
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Average  Error 

Density 

2.076670 

1.67% 

2.094967 

3.03% 

2.172882 

5.94% 

2.287890 

9.25% 

2.539107 

16.01% 

2.855300 

3.090271 

3.227844 

33.87% 

3.329092 

37.90% 

3.427297 

40.97% 

3.814695 

48.93% 

4.202163 

55.16% 

4.533337 

60.26% 

5.056063 

68.17% 

5.464669 

72.73% 

6.321297 

79.80% 

6.990403 

84.22% 

7.913178 

89.40% 

8.625618 

92.77% 

8.970320 

94.59% 

10.121454 

100% 

Figure  43:  Uncertainty  Estimation  in  Yosemite  Scene 
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The  20th  Frame  in  Translating-Tree  Sequence  The  Correct  Flow 

Figure  44:  Translating  Tree  Image  Sequence 

algorithms  at  all  densities. 

Figure  47  shows  the  uncertainty  ellipses  superimposed  on  the  darkened  image.  We 
can  see  that  the  uncertainty  estimation  not  only  detects  the  lack  of  information  along 
certain  direction,  e.g.  the  aperture  problem  along  trunks,  it  also  represented  large 
errors  caused  by  blank  areas  between  trunks. 

A  more  tricky  synthetic  sequence  is  the  diverging  tree  sequence  from  [3,  7]  as 
illustrated  in  Figure  48.  The  significant  difference  between  this  sequence  and  previous 
sequences  is  that  the  underlying  optical  flow  has  a  very  high  gradient,  i.e.  the  image 
velocity  changes  rapidly  as  the  location  in  the  image  changes.  In  fact,  if  we  shift  the 
correct  optical  flow  by  only  two  pixels  horizontally,  the  average  angle  between  the 
shifted  and  the  original  optical  flow  is  1.37°  comparing  to  0.08°  for  the  translating 
tree  sequence  and  0.65°  for  the  Yosemite  sequence.  If  the  window  size  is  larger  than 
two  pixel  widths,  we  can  expect  that  the  average  error  of  the  computed  optical  flow 
should  be  much  larger  than  1.37°.  While  some  results  reported  in  [3]  is  very  close 
or  even  smaller  than  1.37°,  we  believe  that  the  only  reason  those  algorithms  can 
reach  such  a  precision  is  due  to  the  averaging  in  the  temporal  direction.  Actually, 
as  the  authors  pointed  out  in  [3],  the  errors  of  the  differential  approaches  went  up 
substantially  when  only  two  images  are  used.  And  because  the  temporal  averaging  is 
crucial  in  this  sequence,  we  argue  that  the  comparisons  made  in  [3]  (Table  5)  is  not 
fair  for  those  algorithms  with  shorter  temporal  support. 

Using  the  same  parameters  as  those  in  the  translating  tree  sequence,  we  computed 
the  optical  flows  from  adjacent  frames.  Table  8  shows  the  average  errors  when  us¬ 
ing  only  the  20th  and  21st  frame,  and  the  average  errors  when  we  average  fourteen 
computed  optical  flow  fields  from  the  13th  frame  to  the  27th  frame.  First  of  all,  they 
did  verify  our  observation  that  the  average  angular  error  close  to  1.37°  can  only  be 
achieved  by  averaging  along  the  temporal  direction.  Secondly,  the  average  error  of 
our  algorithm  after  averaging  the  fourteen  optical  flow  fields  is  smaller  than  all  the 

results^ll  results  reported  in  [3]  use  at  least  fifteen  frames  except  those  from  Fleet  & 
__ 
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Figure  46:  Average  Error  For  the  Translating  Tree  Scene 
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Average  Angular  Error(°) 

Density 

0.149777 

1.99% 

0.171012 

10.05% 

0.171912 

20.89% 

0.174841 

26.22% 

0.185415 

36.59% 

0.199455 

45.09% 

0.206545 

52.08% 

0.213830 

57.48% 

0.218359 

61.58% 

0.222904 

64.73% 

0.237176 

78.72% 

0.242210 

84.52% 

0.248934 

86.27% 

0.268584 

90.39% 

0.323182 

95.25% 

0.618968 

100% 

Table  7:  Density  vs.  Average  Error  of  Translating  Tree  Scene 
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Figure  47:  Uncertainty  Estimation  in  Translating  Tree  Scene 
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The  20th  Frame  in  Diverging- Tree  Sequence 
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The  Correct  Flow 


Figure  48:  The  Diverging  Tree  Image  Sequence 


Density 

Average  Error  From 
Two  Frames 

Average  Error  From 
Fifteen  Frames 

10% 

2.191357 

1.629412 

20% 

2.188124 

1.595904 

30% 

2.292716 

1.626959 

40% 

2.386003 

1.647144 

50% 

2.483744 

1.696696 

60% 

2.570150 

1.724406 

70% 

2.686011 

1.784787 

80% 

2.825161 

1.863550 

90% 

3.038625 

1.996772 

100% 

3.408894 

2.195650 

Table  8:  Average  Error  of  The  Diverging  Tree  Sequence 


Jepson’s  algorithm  which  achieved  errors  less  than  1.37°.  We  suspect  two  reasons  of 
the  outstanding  performance  of  Fleet  &  Jepson’s  algorithm  reported  in  [3].  One  is 
that  the  extension  of  2D  image  to  3D  spatiotemporal  signal  is  a  better  way  of  tem¬ 
poral  averaging  than  the  sample  averaging  of  the  optical  flows.  The  other  reason  is 
that  the  implementation  of  Fleet  &  Jepson’s  algorithm  utilizes  only  high  frequency 
information  which  is  more  robust  against  large  velocity  gradient. 

We  also  tested  our  algorithm  on  the  four  real  image  sequences  in  [3].  For  all  the 
real  image  sequences,  we  take  two  frames  in  the  middle  of  each  sequence  to  compute 
optical  flows.  The  size  of  the  H  filters  a  =  4.5,  and  the  sampling  limit  in  Figure  39 
is  Ml  =  A^i  =  0,  M2  =  N2  =  15,  which  correspond  to  the  area  in  the  Fourier  domain 
(0.0  <  /r  <  0.86,  —0.86  <  fy  <  0.86).  For  each  experiment,  we  show  four  images: 
the  middle  frame  in  the  sequence,  the  unthresholded  optical  flow,  thresholded  optical 
flow  and  the  uncertainty  estimation. 

The  first  sequence  is  a  scene  of  a  intersection.  Three  cars  are  moving.  The  motions 
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Thresholded  Optical  Flow  Tg  =  0.005  The  Uncertainty  Estimation 

Figure  49:  Optical  Flow  of  The  Taxi  Sequence 

of  the  two  black  cars  are  difficult  to  detect  using  traditional  approaches  since  the 
contrasts  are  low.  The  white  lines  in  the  lower  part  of  the  image  cause  aperture 
problems  which  are  also  detected  in  the  uncertainty  estimation  by  extremely  narrow 
and  long  ellipses.  The  second  sequence  is  a  static  scene  and  the  camera  moves  toward 
the  coke  can  in  the  middle  of  the  image.  The  aperture  problem  and  weak  texture 
problem  happen  almost  everywhere  as  we  can  see  from  the  uncertainty  estimation. 
The  third  sequence  is  a  scene  of  trees  and  the  camera  moves  horizontally.  And  the 
fourth  sequence  is  a  scene  of  a  rotating  Rubic  cube  and  a  platform. 

Though  there  is  no  objective  measure  to  compare  our  results  with  those  reported  in 
[3] ,  we  can  still  conclude  that  in  average  the  performance  of  our  algorithm  is  superior 
to  those  in  [3]  in  both  density  and  precision  by  looking  at  them,  despite  the  fact  that 
we  only  used  two  frames  in  each  sequence. 
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Thresholded  Optical  Flow  =  0.001  The  Uncertainty  Estimation 

Figure  50:  Optical  Flow  of  The  NASA  Coke  Sequence 


The  Seventh  Frame 


Unthresholded  Optical  Flow 
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Thresholded  Optical  Flow  Tg  =  0.001  The  Uncertainty  Estimation 


Figure  51:  Optical  Flow  of  The  SRI  Tree  Sequence 
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Thresholded  Optical  Flow  Tg  =  0.001  The  Uncertainty  Estimation 


Figure  52:  Optical  Flow  of  The  Ruble  Sequence 
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A  Orthogonality  and  Completeness  of  Moment 
Filters 

At  first,  let’s  rephrase  the  theorems  from  [20]  (page  23,  page  29,  page  59)  and  list 
them  as  following: 


Theorem  A.l  The  generalized  Laguerre’s  polynomial  L{n,  a,  x)  is  the  solution  of  the 
following  differential  equation: 

xy"  -1-  (1  -|-  a  —  x)y'  nj/  =  0.  (132) 


Theorem  A. 2  If 

x'=+^+“e-1,=,,6  =  0,(A:  =  0,l,...),  (133) 

at  the  endpoints  of  an  interval  {a,b),  then  the  polynomial  L{n,a,x)  corresponding  to 
different  n’s  are  orthogonal  on  (a, 6),  i.e. 


rb 

/  L{n,  a,  x)L{m,  a,  x)x°‘e~^dx  =  0,  when{m  ^  n). 
J  a 


(134) 


Theorem  A. 3  Let  f{x)  be  continuous  on  a  <  x  <  b  and  have  a  piecewise  continuous 
derivative  in  this  interval,  If  the  integrals 

f  p(x)x‘'e-’dx,  !i[r{x)]^x'*«e-dx 
J  a 

converge,  then  the  following  expansion  is  valid: 

OO 

/(®)  =  ^CnL{n,a,x), 
n—O 

where, 

f!f(^)Hn,a,x)x^e-=^dx 

Cfl  =  — L - - - 

L{n,  a,  x)L{n,  a,  x)x°‘e~^dx 

Let’s  set  /o  =  0  in  Eq.  8,  we  will  prove  that  in  this  case  the  set  of  kfx)  is  orthogonal 
and  complete.  The  case  when  fo^O  can  be  similarly  proved. 

Because  k2n  is  an  even  function,  and  k2m+i  is  an  odd  function,  we  have 


(135) 

(136) 

(137) 


^  (^)  ^  ^2n(^)  ^  0*  (138) 
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And 


2™'^”m!n!  /■+°° 


r°°  rr  ^  i-  X\ 


1  x'^ 


V^cr2(™+”)+i 

2"^+"m!n!  <r  f+°°  1  xr/  1  x  _i  , 

— = — ; - ^ - ;=  /  L(n,~-,t)L{m,—-,t)t  ^dt 

^a^(m+n)+i  ^  Jo  ^  ’  2’  ^  ^  2’  ^ 


da; 


(139) 


where  we  replace  ^  by  i.  Set  the  interval  (a,  b)  to  (0,  +oo),  and  apply  Theorem  A. 2 
to  the  above  equation,  we  have 


Similarly 


<  hnix),  k2m{x)  >—  0,  when(m  ^  n). 


^  ^2n+l  (^)i  ^2m+l  (^)  ^ 

2'"'^"m!n!  f+°°  2r/  1  \  Tr  1  \  -  — 

■“v^<rA™+-)+5  j_^ 

2"i+n™|  I  ^+oo  1  1  j 

m.n.  .^L{m,-,t)th-^dt 


(140) 


dx 


y^^2(m+n)+2  jq 

=  0,  when(m  n).  (141) 

From  Eq.  138,  Eq.  140  and  Eq.  141,  we  prove  the  orthogonality  of  the  moment  filters. 


i.e. 

<  ki{x),kj{x)  >=  0,when(f  ^  j).  (142) 

For  any  function  i^(a;),  it  can  always  be  decomposed  into  sum  of  an  odd  function 
^o(a;)  and  an  even  function  ^e(a;)  as, 

6(^)  =  ^(^W-'^C-^)) 

^e(a;)  =  ^(^(a;)  +  ^(-^)) 

^(x)  =  ^e(a;)  +  6(x).  (143) 

From  Theorem  A. 3,  in  the  interval  (0,+oo),  for  any  function  g{t)  which  satisfies 
Eq.  135,  we  have 


Therefore,  let  t  = 


OO 

—  ^nT(n,  — ,  t), 

n=0  ^ 

(144) 

9{t) 

OO 

=  dmL{m,--,t). 

m-O  ^ 

(145) 

we  have 

= 

OO  1  X^ 

\/2(7  ^  c„xT(n,  2,  ), 

n=0 

(146) 

q  h 

7o|  ^ 

11 

5.2^2)' 

m=0 

(147) 
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Because  g{t)  can  be  any  function,  we  can  let  g{-^)  —  (oi^)  in  Eq.  146,  and  g{-^)  = 
(^e(a^)  in  Eq.  147.  Replacing  them  into  Eq.  143,  in  (0,  oo)  we  have 


i{x) 


+  0O 


1  2  +00  1  2 
y/2a  Y,  Cr,xL{n,  -,  — )  +  Y  dmL{m,  — ) 


n=0 


+  00 

=  E 


Wi 


<  ki{x),ki(x)  > 


2’ 2(72' 
^2(^)5 


m=0 


where,  when  i  =  2m, 


r2m 


W2m  = 


2'"m! 


=  J-r 

j- 


+°°  2'”m! 


1 


1  x'^ 


00  (7 

+00 


2m 


_1 


2’  20-2 


2"^m! 

^/2 

(j2m 

y/irCT  J 

2™m! 

1 

fj2m 

V^<7 

2™m! 

1 

^2m 

\f^a 

<  ((a; 

))  ^2m(' 

2’  2^2  A2(t2  . 

1  X^ 


v2  r°°  ^  i  a;"  , 

-^Jo  ‘-‘dx 


Similarly  when  i  =  2n  +  1 , 

n^2n+i  =<  ^(a;),^2n+i(a;)  > 
Similarly  in  interval  (— oo,0),  we  also  have 

=<  ({x),ki{x)  >  . 


Wi 


(148) 


(149) 


(150) 


(151) 


B  Properties  of  Moment  Filters 

B.l  Recursive  Properties  of  Moment  Filters 

From  Formula  22.3.9  in  [1],  we  have  the  the  explicit  expression  for  the  generalized 
Laguerre  polynomial  as 


n 


L{n,a,x)  =  Y 

771=0 


r(n  +  a+l)(-ir 

r(m  +  Q  +  l)m!(n  —  m)\ 


(152) 
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Therefore,  we  have 

L{n,  1/2,  x)  —  L{n  —  1, 1/2,  x) 


A  r(n  +  3/2)(-l)"  r(n  +  l/2)(-l)-“ 

r(m  +  3/2)m!(n  -  ra)!  ^  r(m  +  3/2)m!(n  -  1  -  m)! 

-A  r(n  +  l/2)(-ir 
^  ^  r(n-l/2  +  l)(-l)’“ 

r(™  -  1/2  +  l)m!(n  -  m)! 

=  T(n,— l/2,x). 

Similarly, 

(n  +  l/2)iy(n,  —1/2,  x)  —  (n  +  l)Z/(n  +  1,  — 1/2,  x)  =  xT(n,  1/2,  x). 
Applying  Eq.  153  and  Eq.  154  to  pi{x)  in  Eq.  2,  we  have 
2"n!  1  x^ 

2”n!  /  1  xy  ,,  1  x’  A 

=  -  2nF2n-i(x)) 


(153) 


(154) 


(155) 


2<r  (  (n  +  l/2)i(n, — )  -  (n  +  l)i(n  +  1,  ,^3) 


2  2o-2- 


=  i((2n  +  l)p2n(x)  -  (TV2n+2(2:)).  (156) 

Actually  Eq.  155  and  Eq.  156  can  be  merged  into  a  single  recursive  formula: 

xpz{x)  =  j{ipi-i{x)  -  cr^p,+i(x)).  (157) 

Replacing  Eq.  157  into  Eq.  1,  we  then  obtain  the  recursive  property  of  the  moment 
filters  in  the  spatial  domain, 

xmi{x)  =  j{imi-i[x)  —  a'^mij^\{x)).  (158) 

Another  recursive  property  of  the  moment  filters  is  with  respect  to  differentiation. 


m'(x)  =  j 


—  j-i _  I  „-jfox 


=  j(mn.i(x)  -  ft,mi(x)). 


(159) 


The  following  recursive  property  of  the  moment  filters  in  the  frequency  domain  is 
obviously  true  from  Eq.  4, 


(/-/o)M.-(/)  =  Mi+,(/). 


(160) 
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B.2  Relations  to  Instantaneous  Frequency  and  Stability  Cri¬ 
terion 

When  we  consider  the  filtering  of  the  whole  signal  instead  of  only  a  specific  spatial 
location,  the  weights  wi  in  Eq.  13  becomes  a  function  of  the  location,  i.e. 

/+00 

({x)mi{x  -  xo)dx.  (161) 

-CO 

We  knew  that  r(;o(a^o)  is  the  output  of  the  Gabor  filter,  let 

u;o(a;o)  ==  Roi^o)  +  jIo{xo), 

where  Ro(xo)  and  /o(a;o)  are  the  real  and  imaginary  parts  of  ryo(a;o)  respectively.  Then 
the  instantaneous  frequency  is  defined  as  the  spatial  derivative  of  phase  value,  i.e. 

RqIq  —  R'qIq 


UXq 


Ri  +  n 

Im(wo«^o) 


Wo 


(162) 


where  Im  means  the  imaginary  part. 

From  Eq.  161  and  Eq.  17,  we  can  compute  the  spatial  derivative  of  wq  as, 

/+00 

^[x)m'o{x  —  xo)dx 

■OO 


Wn 


-OO 
f'+OO 


/±oo 

^{x){Tni{x  -  a:o)  -  fomo{x  -  Xo))dx 

■OO 

=  jifoWo-Wi). 

Replacing  Eq.  163  into  Eq.  162,  we  then  have 

d  r  Re(wiO  ^  „  fwi 

=  fo-  ^  =  /o  -  Re  — 

dxo  Wo  r  \wo 


(163) 

(164) 


where  Re  means  the  real  part. 

The  stability  criterion  of  the  Gabor  filter  output  is  usually  formulated  as  a  thresh¬ 
old  [6,  30], 


where  a  is  the  amplitude. 
Because 


II  1  da 

a  dxo 

1  da 
a  dxo 


+ 


d(j) 

dxo 


—  fo 


(165) 


RqR'o  -h  IqIq 

Rl  +  Ii 

Re(w>o) 

II  ^0  P 

Im(r(;iu;o) 


„  '^0  ,, 

=  -lm(^), 
\WoJ 


(166) 
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we  then  have 


(167) 


II  u;,  IP 
I  ^0  IP’ 


C  Error  Estimations  in  Moment  Filter  Approach 

C.l  Focus 

Suppose  the  energy  of  the  white  noise  in  the  image  is  7V^,  the  noise  energy  in  Ui  or 
Vi  will  be: 

Wf  =  EfNl  (168) 

where  Ef  is  the  energy  of  the  filter  in  Eq.  31. 

Let 


Coin)  =  Uo  +  fonU^-^uil-f^u)U2, 

C,{u)  =  Ur  +  fouU2-~u{l-f^u)U3. 

Replacing  them  back  to  Eq.  39  and  Eq.  40,  and  taking  the  two  error  sources,  i.e. 
white  noise  WQ  and  truncation  error  T(),  into  consideration,  we  have: 

u  =  ^  (ln(Co(u)  +  W{Co)  +  T((7o))  -  \n{Vo  +  vPo  +  W{Vo))) , 

Jo 

+  vPi  +  W{V)  _  C'i(u)  +  lE(C'i)  +  r(C'i) 

Vo  +  vPo  +  W{Vo) ""  Co(u)  +  W{Co)  +  T{Co) ' 

While  the  magnitude  of  the  white  noise  error  W{Co),W{C\),W{Vq)  and  W{Vi) 
can  be  estimated  from  Eq.  168,  the  magnitude  of  the  truncation  error  (assuming  the 
third  order  is  truncated)  can  be  estimated  as 

/7'2 

II  r(Co)  IP  =  II  +  »«,)  -  Co  IP  II  -  /„V)I7„  IP. 

II  nco  IP  =  II  +  vP,)  -  C.  IP  II  -  /„V)Co  IP. 

where  the  first  term  is  the  residual  in  Eq.  39,  and  the  second  the  estimated  energy  of 
truncated  part. 

From  Eq.  169  and  Eq.  170,  assuming  all  the  error  sources  are  independent,  after 
some  manipulations,  we  have. 


(169) 

(170) 


Jl  _ 

2  Co(u) 

{Vo  +  vPo)C[{u)  -  {V  +  vPr)C'o{u) 


Vq  +vPo 

CiPo  -  CoPi 
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1  n  _ _ 1 _  0  — i—  0 

Co(u)  ^  Vo+vPo  ^  Co(u)  ^ 

Vi vP\  — (Vo  +  U-Po)  —Ci{u)  Co{u)  Vi  +  vPi  — (Vo  +  uPo) 


Rewrite  the  above  equation  in  a  more  concise  form,  we  have 


AW, 


where  A  is  a  2x6  complex  matrix,  and  W  is  the  6x1  error  vector. 


W{Co)  ■ 
W{CP) 
W{Vo) 
w{Vi)  • 

T{Co) 

T{Cr)  J 


(171) 


C.2  Stereo 

Suppose  the  energy  of  the  white  noise  in  the  image  is  N^,  then  the  noise  energy  in 
Ui  or  Vi  will  be; 

=  EfNl,  (172) 

where  Ef  is  the  energy  of  the  ith  moment  filter. 

Let 

Co{D)  =  C/o-iPt/i-ip't/2,  (173) 


Ci(P)  =  Ur-jDU2--D^Us, 


(174) 


and  Cq{D)  and  C[{D)  be  their  derivatives  respectively. 

The  magnitude  of  the  truncation  error  (assuming  the  third  order  is  truncated)  can 
be  estimated  as, 

IlnCo)!^  =  ||e-"»'’(V„+,.F„)-C„||^+^||— C7„||^  (175) 


r(c,) 


-  C\  II"  II  —do  II", 


(176) 


where  the  first  term  is  the  residual  in  Eq.  73,  and  the  second  the  estimated  energy  of 
truncated  part. 

We  then  have  the  error  of  D  and  //  as  linear  functions  of  the  two  error  sources,  i.e. 
white  noise  VE()  and  truncation  error  T{),  as  following, 

dD  ]  r  -f/n  — 

_dfx  \  [  (Vo  +  f,Po)C[{D)  -  (El  +  fiPi)C^{D)  CrPo  -  CoPi  \ 

■  VE(C'o)  ■ 
W{Ci) 

®  ^  W(TV)) 

V+fiPi  -{Vo  +  fiPo)  -Ci{D)  Co{D)  Ei+z^Pi  -{Vo  +  fiPo)\  W{Vr) 

T{Co) 

[  nc,)  J 
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Rewrite  the  above  equation  into  a  more  concise  form,  we  have 


=  AW, 


(177) 


where  A  is  a  2x6  complex  matrix,  and  W  is  the  6x1  error  vector. 

D  Properties  of  Hypergeometric  Filters 

In  this  appendix,  we  will  prove  some  properties  of  the  H  filters,  which  are  fundamen¬ 
tal  to  their  applications.  The  properties  include  the  relation  between  H  filters  and 
moment  filters,  the  recursiveness  in  the  frequency  and  spatial  domains,  and  the  re¬ 
cursiveness  with  respect  to  differential  operation.  First  of  all,  let  us  list  the  properties 
of  the  confluent  hypergeometric  function  $(a,c,  a;),  which  are  useful  to  us.  All  the 
properties  listed  below  are  from  Chapter  13  in  [1]: 


$(a,  c,  x)  =  e®$(c  —  a,  c,  —x), 


(178) 


(1  -f  a  ~  c)^{a,c,x)  —  a$(a  +  l,c,x)  +  (c  —  l)#(a,c  —  l,a;)  =  0,  (179) 

c$(a,  c,  x)  —  c$(a  —  1,  c,  x)  —  x^{a,  c  +  1,  x)  =  0,  (180) 


-L{n,a^x)  —  $(— n,a+ 


(181) 


From  Eq.  84  and  Eq.  86,  using  the  relation  between  the  Laguerre’s  polynomial  and 
the  confluent  hypergeometric  function  in  Eq.  181,  we  have 


1  x^ 

h2i{x) h-2i{x)  =  2a„j$(-f, — )e  20-2 

3 

/i2i+i(a;)  - /i_(2;+i)(x)  =  j2h^x<^{-i,-,—)e  2.72 

i\  .  1  x^ 


(182) 


(183) 


in  which  the  right  hand  sides  are  the  moment  filters  (/o  =  0)  multiplied  by  constants. 

From  the  definition  of  the  H  filters  in  Eq.  77,  Eq.  78,  and  Eq.  79,  we  have  the 
recursive  property  of  the  H  filters  in  the  frequency  domain. 


Cm+1 

fHo(f)  = 

^1 

SH.M)  =  — ^/7-(„+.)(/). 

*-m+l 
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To  prove  the  recursive  property  in  the  spatial  domain,  let  us  first  simplify  our 
notations  by  the  following  replacements: 


t  = 


X 


V2a 


(184) 


7rr(m  +  l/2)cr 

y/2  ^  2  ’  ^  2 ’2 


7rr(m  +  l/2)cr 


(185) 


where 


1  —  m  3  9  ' 


=  +  f  )*(i  +  f  4,-<^)i<i86) 

Then  by  using  Eq.  179  and  Eq.  180  and  the  property  of  the  Gamma  function 

r(a;  +  1)  =  xr(a:),  (187) 

we  have  * 

1  ^A+m,,^,l+rn  1  .2,  /^^/2  +  ,2.  ,2  +  m  3  ^2 


2  '.,^,-e)-V2rC-f^m- 


-  -n 

2  ’  2’  ' 


=  -jv^r( 


+J_r(?±!I! 
^  2 


)(*( 


2  +  ml  9,  ^  ,m  1  9 


:7f 


+ 1 i.  -i^) + ,  v/2r(^)*(i:^.  |  -Xx 


—  2  1(0  V^m+l(f)- 


Replacing  Eq.  188  into  Eq.  185,  we  obtain 
jthm(i^  - 


m" 


ITT’  +  1/2  ,  ,  .  I  I, V  ,  ,  . 

^  /im+l(0  -  J  2^  _1 


Then  we  replace  t  in  Eq.  189  by  Eq.  184.  After  some  manipulations,  we  have 


xhra{x)  =  -jo-  I  \lm  +  ^/fm+i(x)  - 


Mm  —  1/2 


1(^)  1  • 


(188) 


(189) 


(190) 
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Similarly 


xho(x)  =  ~ia(hi(x)-h-i{x)), 


(191) 


m(^)  —  3^  y  2 


■/i_(„_i)(x)  .  (192) 


^  m  —  1/2 

About  the  recursive  property  of  H  filters  about  the  differential  operation,  we  have, 

h'Ax)  =  J^-^[nhL{x)]] 


Similarly 


ho{x) 

h'-M 


=  J- 

^'m+l 
^m+1 


j—{hi{x)  -  h_i{x)), 

Cl 

-i— 

^m+1 


(193) 

(194) 

(195) 
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